8
$\begingroup$

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.

$\endgroup$
5
  • $\begingroup$ Related: See stats.stackexchange.com/q/630463/5536 for the simpler $m=1$ case of this question. $\endgroup$
    – a06e
    Commented Jan 26 at 16:03
  • $\begingroup$ By "$\cosh^m(ax)$", do you mean "$(\cosh(ax))^m$" or that $\cosh$ is applied $m$ times? $\endgroup$
    – mhdadk
    Commented Jan 26 at 16:14
  • $\begingroup$ @mhdadk I mean $(\cosh(ax))^m$. $\endgroup$
    – a06e
    Commented Jan 26 at 16:31
  • $\begingroup$ Can you provide values of $a$ and $m$ that would be of interest for you? $\endgroup$
    – Matt F.
    Commented Jan 27 at 13:42
  • $\begingroup$ $m$ is a positive integer, and $a$ a positive real. @MattF. $\endgroup$
    – a06e
    Commented Jan 27 at 16:46

3 Answers 3

11
$\begingroup$

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.

$\endgroup$
10
  • 3
    $\begingroup$ This is fine for integral $m$ -- but that's the easy case. $\endgroup$
    – whuber
    Commented Jan 26 at 19:56
  • 2
    $\begingroup$ @a06e For integer $m$, the mean of each component is $-a(m-2k)$. $\endgroup$ Commented Jan 27 at 17:08
  • 1
    $\begingroup$ This, then, is a good answer +1. $\endgroup$
    – whuber
    Commented Jan 27 at 18:55
  • 1
    $\begingroup$ +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}$$ $\endgroup$
    – JimB
    Commented Jan 27 at 19:55
  • 1
    $\begingroup$ The $\exp(\mu_k^2)$ in the last equation above should of course be $\exp(\mu_k^2/2)$. $\endgroup$
    – JimB
    Commented Jan 27 at 22:08
6
$\begingroup$

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

$\endgroup$
2
  • $\begingroup$ Ah I almost missed this wonderful answer, thanks! $\endgroup$
    – a06e
    Commented Jan 27 at 19:11
  • 2
    $\begingroup$ @a06e Thanks, this was a neat problem! $\endgroup$ Commented Jan 27 at 19:30
3
$\begingroup$

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

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.