Generating random variable which has a power distribution of Box and Tiao (1962)

Question

Box and Tiao (Biometrika 1962) use a distribution whose density has the following form: $$f(x; \mu, \sigma, \alpha) = \omega \exp\left\{ -\frac{1}{2} \Big\vert\frac{x-\mu}{\sigma}\Big\vert^{\frac{2}{(1+\alpha)}} \right\},$$ where $\omega^{-1} = [\Gamma(g(\alpha)]\,2^{g(\alpha)}\sigma$ is the normalizing constant with $g(\alpha) = \frac{3}{2} + \frac{\alpha}{2},$ $\sigma \gt 0,$ and $-1 \lt \alpha \lt 1$.

When $\alpha=0$ this reduces to the normal distribution; when $\alpha=1$ it reduces to the double exponential (Laplace) distribution, and when $\alpha \to -1^{+}$ it tends to a uniform distribution.

How can I generate random numbers from this distribution for any such value of $\alpha$?

Answer 1

Box & Tiao refer to this as a "convenient class of power distributions," referencing Diananda (1949), Box (1953), and Turner (1960).

Because $\mu$ and $\sigma$ just establish a unit of measurement and the absolute value reflects values around the origin, the basic density is proportional to $\exp(-z^p/2)$ where $p = 2/(1+\alpha)$ and $z \ge 0.$ Changing variables to $y = z^p$ for $0\lt p \lt \infty$ changes the probability element to

$$\exp(-z^p/2)\mathrm{d}z \to \exp(-y/2) \mathrm{d}\left(y^{1/p}\right) = \frac{1}{p}y^{1/p - 1}e^{-y/2}\mathrm{d}y.$$

Since $ p = 2/(1+\alpha),$ this is proportional to a scaled Gamma$(1/p)$ = Gamma$((1+\alpha)/2)$ density, also known as a Chi-squared$(1+\alpha)$ density.

Thus, to generate a value from such a distribution, undo all these transformations in reverse order:

Generate a value $Y$ from a Chi-squared$(1+\alpha)$ distribution, raise it to the $2/(1+\alpha)$ power, randomly negate it (with probability $1/2$), multiply by $\sigma,$ and add $\mu.$

This R code exhibits one such implementation. n is the number of independent values to draw.

rf <- function(n, mu, sigma, alpha) {
  y <- rchisq(n, 1 + alpha)                # A chi-squared variate
  u <- sample(c(-1,1), n, replace = TRUE)  # Random sign change
  y^((1 + alpha)/2) * u * sigma + mu
}

Here are some examples of values generated in this fashion (100,000 of each) along with graphs of $f.$

Generating Chi-squared variates with parameter $1+\alpha$ near zero is problematic. You can see this code works for $1+\alpha = 0.1$ (bottom left), but watch out when it gets much smaller than this:

The spike and gap in the middle should not be there.

The problem lies with floating point arithmetic: even double precision does not suffice. By this point, though, the uniform distribution looks like a good approximation.

---

Appendix

This R code produced the plots. It uses the showtext library to access a Google font for the axis numbers and labels. Few of these fonts, if any, support Greek or math characters, so I had to use the default font for the plot titles (using mtext). Otherwise, everything is done with the base R plotting functions hist and curve. Don't be concerned about the relatively large simulation size: the total computation time is far less than one second to generate these 400,000 variates.

library(showtext)
if(!("Informal" %in% font_families())) font_add_google("Fuzzy Bubbles", "Informal")
showtext_auto()
#
# Density calculation.
#
f <- function(x, mu, sigma, alpha) 
              exp(-1/2 * abs((x - mu) / sigma) ^ (2 / (1 + alpha)))
C <- function(mu, sigma, alpha, ...) 
              integrate(\(x) f(x, mu, sigma, alpha), -Inf, Inf, ...)$value
#
# Specify the distributions to plot.
#
Parameters <- list(list(mu = 0, sigma = 1, alpha = 0),
                   list(mu = 10, sigma = 2, alpha = 1/2),
                   list(mu = 0, sigma = 3, alpha = -0.9),
                   list(mu = 0, sigma = 4, alpha = 0.99))
#
# Generate the samples and plot summaries of them.
#
n.sim <- 1e5 # Sample size per plot
set.seed(17) # For reproducibility
pars <- par(mfrow = c(2, 2), mai = c(1/2, 3/4, 3/8, 1/8)) # Shrink the margins
for (parameters in Parameters)
  with(parameters, {
    x <- rf(n.sim, mu, sigma, alpha)
    hist(x, freq = FALSE, breaks = 100, family = "Informal", 
         xlab = "", main = "", col = gray(0.9), border = gray(0.7))
    mtext(bquote(list(mu==.(mu), sigma==.(sigma), alpha==.(alpha))), 
          cex = 1.25, side = 3, line = 0)
    omega <- 1 / C(mu, sigma, alpha) # Compute the normalizing constant
    curve(omega * f(x, mu, sigma, alpha), add = TRUE, lwd = 2, col = "Red")
  })
par(pars)