Proofs from the notebook Olivier's open lab book.

/ / / / /

Two Sampling Algorithms for Trigonometric Densities

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

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

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

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

where

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

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:

  1. Let $(x, y)​$ be uniformly distributed over $[0, 2\pi) \times [0, C]​$;
  2. 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

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

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.

  1. Let $x \sim f_a $.
  2. 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.