8

Is there an efficient algorithm to draw samples $x \sim P(x)$ from the PDF:

$$ P(x) \propto \cosh^{m}(a x) e^{-x^{2}/2} $$

where $a\ge0$ is a real parameter, and $m$ a positive integer?

Since this is a univariate PDF, my first attempt was to compute the CDF and use the inverse CDF method to sample it, but neither the CDF nor the inverse CDF seem to have a tractable analytical form.

Jarle Tufto
  • 10,939
a06e
  • 4,410
  • 1
  • 22
  • 50

3 Answers3

11

The distribution that the OP seeks is a straightforward extension of the formula for the case $m=1$ that Xi'an derived in the linked answer.

Suppose that $P(x) \propto (\cosh(ax))^m \exp(-x^2/2)$. Then, since \begin{align}(\cosh(ax))^m &= \left(\frac{\exp(ax)+\exp(-ax)}{2}\right)^m\\ &= \sum_{k=0}^m 2^{-m} \binom{m}{k}\exp(akx)\exp(-a(m-k)x)\\ &= \sum_{k=0}^m 2^{-m} \binom{m}{k} \exp(-a(m-2k)x), \end{align} $P(x)$ is once again a weighted sum (also called a mixture) of normal distributions with means ranging from $-ma$ to $+ma$ in steps of $2a$. Note that the standard normal distribution is included in the mixture only when $m$ is even. The weights are given by the probability mass distribution function of a $(m,\frac 12)$ binomial random variable.

Turning to the actual question asked by the OP about an efficient algorithm for sampling from this distribution, consider the various answers to this question, some of which even contain R code with varying degrees of efficiency.

Dilip Sarwate
  • 46,658
  • 3
    This is fine for integral $m$ -- but that's the easy case. – whuber Jan 26 '24 at 19:56
  • If I understand correctly, sampling $P(x)$ then amounts to drawing $k \sim \operatorname{Binom}(m,1/2)$ (which is easy), and then $x \sim \mathcal N(2ka, 1)$ (which is also easy) ? – a06e Jan 27 '24 at 11:32
  • @whuber I am interested in integer $m$ – a06e Jan 27 '24 at 16:47
  • 2
    @a06e For integer $m$, the mean of each component is $-a(m-2k)$. – Jarle Tufto Jan 27 '24 at 17:08
  • 1
    This, then, is a good answer +1. – whuber Jan 27 '24 at 18:55
  • @JarleTufto Right thanks. – a06e Jan 27 '24 at 19:10
  • 1
    +1 but don't the weights (prior to making them sum to 1) need to depend on $\mu_k=-a(m-2k)$? $$\begin{align}(\cosh(ax))^m \exp(-x^2/2) &= \sum_{k=0}^m 2^{-m} \binom{m}{k} \exp(-a(m-2k)x) \exp(-x^2/2) \ &= \sum_{k=0}^m 2^{-m} \binom{m}{k} \exp(\mu_k x-x^2/2) \ &= \sum_{k=0}^m 2^{-m} \binom{m}{k} \exp(\mu_k x-x^2/2 - \mu_k^2/2 + \mu_k^2/2) \ &= \sum_{k=0}^m 2^{-m} \binom{m}{k} \exp(\mu_k^2)\exp(-(x-\mu_k)^2/2) \end{align}$$ – JimB Jan 27 '24 at 19:55
  • 1
    The $\exp(\mu_k^2)$ in the last equation above should of course be $\exp(\mu_k^2/2)$. – JimB Jan 27 '24 at 22:08
  • 1
    @JimB $P(x)$, the pdf of $X$ is proportional to $\sum_{k=0}^m 2^{-m} \binom{m}{k} \exp(-a(m-2k)x)\exp(-x^2/2$, not equal to $\sum_{k=0}^m 2^{-m} \binom{m}{k} \exp(-a(m-2k)x)\exp(-x^2/2)$. The $\exp(\mu_k^2/2)$ factors that you have discovered are just part of the constant of proportionality. The weights $2^{-m}\binom{m}{k}$ do sum to $1$. – Dilip Sarwate Jan 27 '24 at 23:03
  • 1
    I totally agree with your proportionality comment but $\exp(\mu_k^2/2)$ varies for every value of $k$ and therefore isn't a constant part of the constant of proportionality. The weights (prior to dividing by the total of the weights) does need to include $\exp(\mu_k^2/2)$. With values of $a>0$ and even relatively small values of $m$, the two extreme normal distributions (associated with $k=0$ and $k=m$) completely dominate with all of the other standardized weights nearly zero. – JimB Jan 28 '24 at 19:38
6

Building on the answer of Dilip, for non-integer $m$, letting $\{m\}=m-\lfloor m\rfloor$ denote the fractional part of $m$, and using the fact that $\cosh x\le \exp|x|$, the (unnormalized) target density can be rewritten as $$ \begin{align} f(x)&=\cosh^m(ax)\exp(-x^2/2) \\&=\cosh^{\{m\}}(ax)\cosh^{\lfloor m\rfloor} (ax)\exp(-x^2/2) \\&\le \exp(a\{m\}|x|)\cosh^{\lfloor m\rfloor} (ax)\exp(-x^2/2) \\&= \sum_{k=0}^{\lfloor m\rfloor} 2^{-{\lfloor m\rfloor}} \binom{{\lfloor m\rfloor}}{k} \exp\Big(-x^2/2-a({\lfloor m\rfloor}-2k)x+a\{m\}|x|\Big). \end{align} $$ Like the target $f(x)$, the right hand side is symmetric around zero. For $x\ge 0$, it is proportional to the density of a mixture of $\lfloor m\rfloor+1$ Gaussian distributions truncated below zero. Completing the square in $x$, the mean (before truncation) of component $k$ is $$ \mu_k=-a(\lfloor m\rfloor-\{m\}-2k) $$ and the associated weights are $$ w_k \propto \Phi(\mu_k)\binom{\lfloor m\rfloor}{k}e^{\mu_k^2/2}. $$

We can easily simulate proposals efficiently from this mixture (using the alias and inversion methods), and then obtain samples from $f(x)$ by accepting the proposals with probabilities given by $f(x)$ divided by the above right hand side.

The acceptance probability simplifies to $$ \left(\frac{1+\exp(-2a|x|)}{2}\right)^{\{m\}}>\frac1{2^{\{m\}}} $$ so the overall acceptance probabilty is always greater than 1/2.

R implementation:

dcoshgauss <- function(x, a, m)  
  cosh(a*x)^m*exp(-x^2/2)
rcoshgauss <- function(n, a, m) {
  floorm <- floor(m)
  fractm <- m - floorm
  kk <- 0:floorm
  mu <- -a*(floorm - fractm - 2*kk)
  w <- pnorm(mu)*choose(floorm, kk)*exp(mu^2/2)
  w <- w/sum(w)
  xx <- numeric(n)
  i <- 0
  s <- 0
  while (i < n) {
    k <- sample(kk+1, size = 1, prob=w)
    u <- runif(1, pnorm(-mu[k]), 1)
    x <- mu[k] + qnorm(u)
    acceptprob <- ((1+exp(-2*x*a))/2)^fractm
    u <- runif(1)
    s <- s + 1
    if (u<acceptprob) {
      sign <- sample(c(-1,1), 1)
      x <- sign*x
      i <- i + 1
      xx[i] <- x
    }
  }
  cat("Acceptance probability was",n/s,"\n")
  xx
}

set.seed(1)
x <- rcoshgauss(1e6, 1, 2.9)
#> Acceptance probability was 0.5489607
hist(x, breaks=200, xlim=c(-8,8), main="")


curve(dcoshgauss(x, 1, 2.9), -8,8)

Created on 2024-01-28 with reprex v2.0.2

Jarle Tufto
  • 10,939
2

This is just an extended comment.

With the restriction of the values of $m$ being positive integers there is an analytic solution of a mixture of normals each with variance 1 as shown by @DilipSarwate (although I currently disagree about the weights). However, showing graphic examples should almost always be considered in which interesting features can appear.

In this case when $a$ and $m$ are large enough (and they don't need to be very large) one is essentially sampling from a mixture of just two normals with means of opposite sign and variance 1. This might also be the case for large values of $m$ when $m$ is not an integer.

Below are a few examples showing the resulting density (blue) and the individual contributions (red).

Small $a$ and small $m$ Small a and small m

$a=1$ and $m=2$ a=1 and m=2

$a=1$ and $m=3$ a=1 and m=3

$a=1$ and $m=4$ a=1 and m=4

JimB
  • 3,734
  • 11
  • 20