Trigonometric densities (or non-negative trigonometric sums) are probability density functions of circular random variables (i.e. \(2\pi\)-periodic densities) which take the form \[ f(u) = a_0 + \sum_{k=1}^n(a_k \sin(k u) + b_k\cos(ku)) \tag{1} \]
for some real coefficients \(a_k, b_k \in \mathbb{R}\) which are such that \(f(u) \geq 0\) and \(a_0 = \frac{1}{2\pi} \int f(u)\,du = (2\pi)^{-1}\). These provide flexible models of circular distributions. Circular density modelling comes up in studies about the mechanisms of animal orientation and also come up in bio-informatics in relationship to the protein structure prediction problem (the secondary structure of a protein - the way its backbone folds - is determined by a sequence of angles).
Here I am discussing two simple sampling algorithms for such trigonometric densities. The first is the rejection sampling algorithm proposed in Fernández-Durán et al. (2014) and the second uses negative mixture sampling.
Parametrizing trigonometric densities
By Féjer’s Theorem, the conditions on the coefficients \(a_k\) and \(b_k\) can be stated as follows: there exists a vector of complex coefficients \(c = (c_0, c_1, \dots, c_n)\) with \(\|c\|^2 = (2\pi)^{-1}\) and satisfying
\[ f(u) = \left\| \sum_{k=0}^n c_k e^{ik u} \right\|^2. \tag{2} \]
This provides an explicit parametrization of the space of trigonometric densities in terms of a complex hypersphere. See Fernandez-Duran (2004) for more details.
Density basis of the trigonometric polynomials
In Binette & Guillotte (2019), we studied the De la Vallée Poussin density basis of the trigonometric polynomials given by
\[ C_{j,n}(u) = \frac{2^n}{2\pi {2n \choose n}} \left(1+\cos\left(u - \tfrac{2\pi j}{2n+1}\right)\right)^n,\quad j\in \{0,1,\dots, 2n\}. \tag{3} \]
These can be used to express trigonometric densities as mixtures of probability density functions (instead of the functions \(\cos\) and \(\sin\), and the change of basis formula follows from the expression
\[ C_{j,n}(u) = T_{j,n}\,\left[e^{-i nu}\; \cdots\; e^{-i u}\; 1\; e^{i u}\; \cdots\; e^{i nu}\right]^{T} \]
where
\[ T_{j,n} = \left[\exp\left\{ -i\frac{2\pi j p}{2n+1} u \right\}{2n \choose n-p} \Big / {2n \choose n}\right]_{p \in \{-n, \dots, n\}}. \]
We’re using the complex functions \(e^{i2\pi k u}\) instead of \(\sin\) and \(\cos\) simply because they are neater to work with; it doesn’t change much otherwise.
We also show in our paper that if \(V \sim \text{Ber}(1 / 2)\) and \(W \sim \text{Beta}(1 / 2, 1 / 2+n)\), then
\[ (1-2V)\arccos(1-2W) +\tfrac{2\pi j}{2n+1} \sim C_{j,n}. \]
This provide an easily formula to sample from the basis functions \(C_{j,n}\) and their mixtures.
Algorithm 1: Naive rejection sampling
Given an uniform upper bound \(C\) on the family \(\mathcal{V}_n\) of trigonometric densities, we can sample from a given \(f\in \mathcal{V}_n\) using simple rejection sampling as follows:
- Let \((x, y)\) be uniformly distributed over \([0, 2\pi) \times [0, C]\);
- If \(y \leq f(x)\), then return \(x\); otherwise return to step 1.
Now the problem is to figure out a good upper bound \(C\). The most basic idea is to do as in Fernandez-Duran et al. (2014) and to apply the Cauchy-Schwarz inequality
\[ f(u) = \left\| \sum_{k=0}^n c_k e^{i k u} \right\|^2 \leq \|c\|^2 \sum_{k=0}^n|e^{iku}| = \frac{n+1}{2\pi}. \]
Can we find a better bound? I think that \(C = \sqrt{n}/\pi\) would work, but I have no clue how to prove it….
Let’s implement this in R.
Implementation
First we need a trigonometric density model.
trig_function <- function(c_real, complex=NULL) {
# Returns the trigonometric function defined as either:
# f(u) = 1/(2\pi) + \sum_{k=1}^{n} c_real[2*k-1] \sin(k u) + c_real[2*k] \cos(ku),
# or
# f(u) = \| \sum_{k=0}^n complex e^{i k u} \|^2,
# where n is the degree of the polynomial.
#
# Args
# c_real: Vector of 2*n real numbers, where n is the degree of
# the trigonometric polynomial.
# complex: Vector of (n+1) complex numbers.
if (!is.null(complex)) {
lambd <- function(u) {
n = length(complex) - 1
k = 0:n
return(abs(sum(complex * exp(u * k * 1i)))**2)
}
}
else {
lambd <- function(u) {
n = length(c_real)/2
k = 1:n
return(1/(2*pi) + sum(c_real[2*k - 1] * cos(k*u)) + sum(c_real[2*k] * cos(k*u)))
}
}
return(Vectorize(lambd));
}We can also generate random trigonometric densities of a fixed degree as follows.
rtrig <- function(n) {
u = rnorm(n);
v = rnorm(n);
c_comp = u + v*1i;
c_comp = c_comp / (sqrt(2*pi*sum(abs(c_comp)**2)));
return(trig_function(complex=c_comp))
}Usage is like this:
u = seq(0, 2*pi, 0.005)
plot(u, rtrig(10)(u), type="l")
And finally we can implement the naive rejection sampling algorithm.
naive_rejection_sampling <- function(f, n) {
# Returns a random variate following the trigonometric density f of degree n.
drawn = FALSE
while(!drawn) {
x = runif(1)*2*pi
y = runif(1)*(n+1) / (2*pi)
if (y < f(x)) {
drawn = TRUE
}
}
return(x);
}Algorithm 2: Negative Mixture Sampling
Another approach to simulate from trigonometric densities relies on the De la Vallée Poussin mixture representation. That is, any \(f\in \mathcal{V}_n \) can be written as
\[ f = \alpha f_a - (\alpha - 1) f_b,\qquad f_a = \sum_{j=0}^{2n} a_j C_{j_n}, \quad f_b = \sum_{j=0}^{2n} b_j C_{j,n}, \]
where \(\alpha \geq 1 \), \(a_j, b_j \geq 0 \) and \(\sum_ j a_j = \sum_ j b_ j = 1\). We can assume that \(a_j b_j = 0 \) for every \(j \); i.e. there is no redundancy in the components of \(f_a \) and \(f_b \). The density \(f_b \) accounts for negative weights in the mixture representation of \(f \) using the De la Vallée Poussin densities \((3) \).
We can now sample from \(f \) using samples from \(f_a \) and a simple rejection method.
Algorithm 2.
- Let \(x \sim f_a \).
- Return \(x \) with probability \(\frac{f(x)}{\alpha f_a(x)} \); otherwise return to step 1.
Implementation
De la Vallée Poussin densities and its random variate generator.
dvallee <- function(u, j, n) {
# De la Vallée Poussin density $C_{j,n}(u)$
return(2^n * (1+cos(u - (2*pi*j)/(2*n+1)))^n / (2*pi*choose(2*n, n)))
}rvallee <- function(j, n, m) {
# Returns m random variates following the De la Vallée Poussin density $C_{j,n}$.
V = runif(m) > 0.5
W = rbeta(m, 1/2, 1/2 + n)
return((1-2*V)*acos(1-2*W) + (2*pi*j)/(2*n + 1))
}Usage:
s = rvallee(2, 5, 10000)
u = seq(-pi, pi, 0.05)
hist(s, prob=TRUE, xlim=c(-pi,pi))
lines(u, dvallee(u, 2, 5), col=2)
De la Vallée Poussin mixtures.
dValleeMixture <- function(coeffs) {
# De la Vallée Poussin mixture densities
n = (length(coeffs) - 1)/2;
lambd <- function(u) {
j = 0:(2*n)
return(sum(dvallee(u, j, n) * coeffs))
}
return(Vectorize(lambd))
}rValleeMixture <- function(coeffs) {
# Random sample from a De la Vallée Poussin mixture density. The mixture weights are allowed to take negative values.
f = dValleeMixture(coeffs)
n = (length(coeffs) - 1)/2
a = coeffs * (coeffs > 0)
b = coeffs * (coeffs < 0)
alpha = sum(a)
a = a / alpha
b = b / (1-alpha)
fa = dValleeMixture(a)
drawn = FALSE
while(!drawn) {
# Sample from f_a
i = sample(0:(2*n), 1, prob = a)
x = rvallee(i, n, 1)
if ( runif(1) < f(x)/(alpha*fa(x))) {
drawn = TRUE
}
}
return(x %% (2*pi))
}Example:
coeffs = c(0.55, -0.15, 0.55, 0, 0, 0,0.05)
f = dValleeMixture(coeffs)
u = seq(0, 2*pi, 0.05)
s = replicate(50000, rValleeMixture(coeffs))
hist(s, prob=T, ylim=c(0, 0.6))
lines(u, f(u), col=2)
Other things we could do:
- The black box Lipschitz sampling algorithm can also be used to sample from trigonometric densities. This requires to compute good upper bounds on the Lipchitz constant on the density, which should be doable using the De la Vallée Poussin mixture representation.