8
$\begingroup$

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

$\endgroup$
0

1 Answer 1

10
$\begingroup$

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

enter image description here

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: enter image description here

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)
$\endgroup$
7
  • 3
    $\begingroup$ That's some clean code... $\endgroup$
    – Zen
    Oct 5, 2022 at 17:20
  • 1
    $\begingroup$ Nice and very fastly delivered:) answer. $\endgroup$
    – Yves
    Oct 5, 2022 at 17:32
  • $\begingroup$ Beautiful solution. Thank you! $\endgroup$
    – user67724
    Oct 6, 2022 at 18:19
  • $\begingroup$ @whuber: Can you please show us how you generated the lovely plots? $\endgroup$
    – user67724
    Oct 6, 2022 at 18:26
  • $\begingroup$ @user67724 Done. $\endgroup$
    – whuber
    Oct 6, 2022 at 19:05

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.