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

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$$?

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)
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)

• That's some clean code...
– Zen
Oct 5, 2022 at 17:20
• Nice and very fastly delivered:) answer.
– Yves
Oct 5, 2022 at 17:32
• Beautiful solution. Thank you! Oct 6, 2022 at 18:19
• @whuber: Can you please show us how you generated the lovely plots? Oct 6, 2022 at 18:26
• @user67724 Done.
– whuber
Oct 6, 2022 at 19:05