Working with the Exponential Power Distribution Using gnorm

Exponential Power Density, dgnorm

An exponential power distributed random variable, \(x\), has the following density:

\[ p(x | \mu, \alpha, \beta) = \frac{\beta}{2\alpha \Gamma(1/\beta)}\text{exp}\{-(\frac{|x - \mu|}{\alpha})^\beta\}, \] where \(\mu\) is the mean of \(x\) and \(\alpha > 0\) and \(\beta > 0\) are scale and shape parameters.

The function dgnorm gives the exponential power density given above, \(p(x | \mu, \alpha, \beta)\), evaluated at specified values of \(x\), \(\mu\), \(\alpha\) and \(\beta\). Below, we give an example of evaluating the exponential power density using dgnorm for \(\mu = 0\), \(\alpha = \sqrt{2}\) and \(\beta = 2\). Note that this is in fact the standard normal density!

xs <- seq(-1, 1, length.out = 100)
plot(xs, dgnorm(xs, mu = 0, alpha = sqrt(2), beta = 2), type = "l", 
     xlab = "x", ylab = expression(p(x)))

Like dnorm, dgnorm has a log argument, with the default log=FALSE. When log=TRUE is specified, \(\text{log}p(x | \mu, \alpha, \beta)\) is returned.

Exponential Power CDF, pgnorm

The exponential power CDF is derived as follows: \[ \begin{aligned} \text{Pr}(x \leq q | \mu, \alpha, \beta) &= \int_{-\infty}^q \frac{\beta}{2\alpha \Gamma(1/\beta)}\text{exp}\{-(\frac{|x - \mu|}{\alpha})^\beta\} d x \\ &=\int_{-\infty}^{q - \mu} \frac{\beta}{2\alpha \Gamma(1/\beta)}\text{exp}\{-(\frac{1}{\alpha})^\beta |z|^\beta \} d z \\ &=\Bigg\{\begin{array}{cc} \int_{-\infty}^{q - \mu} \frac{\beta}{2\alpha \Gamma(1/\beta)}\text{exp}\{-(\frac{1}{\alpha})^\beta \left(-z\right)^\beta \} d z & q \leq \mu \\ \int_{-\infty}^{0} \frac{\beta}{2\alpha \Gamma(1/\beta)}\text{exp}\{-(\frac{1}{\alpha})^\beta (-z)^\beta + \int_{0}^{q - \mu} \frac{\beta}{2\alpha \Gamma(1/\beta)}\text{exp}\{-(\frac{1}{\alpha})^\beta z^\beta \} d z & q > \mu \\ \end{array} \\ &=\Bigg\{\begin{array}{cc} \int_{-(q - \mu)}^{\infty} \frac{\beta}{2\alpha \Gamma(1/\beta)}\text{exp}\{-(\frac{1}{\alpha})^\beta z^\beta \} d z & q \leq \mu \\ \int_{0}^{\infty} \frac{\beta}{2\alpha \Gamma(1/\beta)}\text{exp}\{-(\frac{1}{\alpha})^\beta z^\beta dz + \int_{0}^{q - \mu} \frac{\beta}{2\alpha \Gamma(1/\beta)}\text{exp}\{-(\frac{1}{\alpha})^\beta z^\beta \} dz & q > \mu \\ \end{array} \\ &=\Bigg\{\begin{array}{cc} \int_{0}^{\infty} \frac{\beta}{2\alpha \Gamma(1/\beta)}\text{exp}\{-(\frac{1}{\alpha})^\beta z^\beta \} d z - \int_{0}^{-(q - \mu)} \frac{\beta}{2\alpha \Gamma(1/\beta)}\text{exp}\{-(\frac{1}{\alpha})^\beta z^\beta \} d z& q \leq \mu \\ \int_{0}^{\infty} \frac{\beta}{2\alpha \Gamma(1/\beta)}\text{exp}\{-(\frac{1}{\alpha})^\beta z^\beta\} dz+ \int_{0}^{q - \mu} \frac{\beta}{2\alpha \Gamma(1/\beta)}\text{exp}\{-(\frac{1}{\alpha})^\beta z^\beta \} d z & q > \mu \\ \end{array} \\ &= \int_{0}^{\infty} \frac{\beta}{2\alpha \Gamma(1/\beta)}\text{exp}\{-(\frac{1}{\alpha})^\beta z^\beta \}dz+ \text{sign}(q - \mu) \int_{0}^{|q - \mu|} \frac{\beta}{2\alpha \Gamma(1/\beta)}\text{exp}\{-(\frac{1}{\alpha})^\beta z^\beta \} d z \\ &= \int_{0}^{\infty} \frac{w^{1/\beta - 1}}{2\alpha \Gamma(1/\beta)}\text{exp}\{-(\frac{1}{\alpha})^\beta w \}dw+ \text{sign}(q - \mu) \int_{0}^{|q - \mu|^\beta} \frac{w^{1/\beta - 1}}{2\alpha \Gamma(1/\beta)}\text{exp}\{-(\frac{1}{\alpha})^\beta w \} d w && z^\beta = w \\ &= \frac{1}{2}+ \frac{\text{sign}(q - \mu)}{2} \underbrace{\int_{0}^{|q - \mu|^\beta} \frac{w^{1/\beta - 1}}{\alpha \Gamma(1/\beta)}\text{exp}\{-(\frac{1}{\alpha})^\beta w \} d w}_{\text{Gamma CDF, Shape=$\frac{1}{\beta}$, Rate=$(\frac{1}{\alpha})^\beta$}} \end{aligned} \] We can evaluate the exponential power CDF using the gamma CDF - the pgamma function evaluates the CDF of a gamma distribution with parameters \(\frac{1}{\beta}\) and \((\frac{1}{\alpha})^\beta\) at \(|q - \mu|^\beta\).

The function pgnorm returns the value of the CDF for specified values of \(q\), \(\mu\), \(\alpha\) and \(\beta\). Like the corresponding pnorm function, it also has the following arguments:

• log.p, which returns the log probability when set to TRUE. The default is log.p=FALSE;
• lower.tail, which returns \(\text{Pr}(x > q | \mu, \alpha, \beta)\) when set to FALSE. The default is lower.tail=TRUE.

Below, we give an example of evaluating the exponential power CDF using pgnorm for \(\mu = 0\), \(\alpha = \sqrt{2}\) and \(\beta = 2\). As in the previous example - this is the standard normal CDF.

xs <- seq(-1, 1, length.out = 100)
plot(xs, pgnorm(xs, 0, sqrt(2), 2), type = "l", xlab = "q", ylab = expression(paste("Pr(", x<=q, ")", sep = "")))

Exponential Power Quantile Function (Inverse CDF), qgnorm

A (relatively) straightforward way to compute the inverse CDF function for an exponential power random variable is to use its scale-sign representation (Gupta and Varga, 1993)⁴. We can write \(x \stackrel{d}{=} su + \mu\), where \(s\) and \(u\) are independent random variables with \(p(s | \alpha, \beta) = \frac{\beta}{\alpha \Gamma(1/\beta)}\text{exp}\{-(\frac{s}{\alpha})^\beta\}\) and \(u\) uniformly distributed on \(\{-1, 1\}\).

Given \(p\), we want to find the value of \(q\) that satisfies \(p = \text{Pr}(x \leq q | \mu, \alpha, \beta)\). We can rewrite the problem using \(s\) and \(u\) as finding \(p\) and \(q\) that satisfy: \[ \begin{aligned} |p - 0.5|&= \text{Pr}(s \leq |q - \mu| | \alpha, \beta)\text{Pr}(u =\text{sign}(q - \mu)) \\ &= \frac{1}{2}\text{Pr}(s \leq |q - \mu| | \alpha, \beta) \\ &= \frac{1}{2}\int_0^{|q - \mu|} \frac{\beta}{\alpha \Gamma(1/\beta)}\text{exp}\{-(\frac{s}{\alpha})^\beta\} ds \\ &= \frac{1}{2}\underbrace{\int_0^{|q - \mu |^\beta} \frac{1}{\alpha \Gamma(1/\beta)}t^{1/\beta - 1}\text{exp}\{-(\frac{1}{\alpha})^\beta t\} dt}_{\text{Gamma CDF, Shape=$\frac{1}{\beta}$, Rate=$(\frac{1}{\alpha})^\beta)$}} & s^\beta = t \end{aligned} \]

Let \(g(2|p - 0.5|; \frac{1}{\beta}, (\frac{1}{\alpha})^\beta)\) refer to the inverse CDF function for a gamma distribution with shape \(\frac{1}{\beta}\) and rate \((\frac{1}{\alpha})^\beta\) evaluated at \(2|p - 0.5|\). Then the inverse CDF of the exponential power distribution is given by: \[ |q - \mu| = g(2|p - 0.5|; \frac{1}{\beta}, (\frac{1}{\alpha})^\beta)^{1/\beta}. \]

Noting that \(q > \mu\) when \(p > 0.5\) and \(q \leq \mu\) when \(p \leq 0.5\), we get: \[ q = \text{sign}(p - 0.5)g(2|p - 0.5|; \frac{1}{\beta}, (\frac{1}{\alpha})^\beta)^{1/\beta} + \mu. \]

This is what is implemented by qgnorm for specified values of \(q\), \(\mu\), \(\alpha\) and \(\beta\). Like the corresponding qnorm function, it also has the following arguments:

• log.p, which indicates that the log probability has been provided when set to TRUE. The default is log.p=FALSE;
• lower.tail, which indicates that \(q\) satisfies \(p \text{Pr}(x > q | \mu, \alpha, \beta)\) when set to FALSE. The default is lower.tail=TRUE.

Below, we give an example of evaluating the exponential power inverse CDF using qgnorm for \(\mu = 0\), \(\alpha = \sqrt{2}\) and \(\beta = 2\). As in the previous examples - this is the standard normal CDF.

xs <- seq(0, 1, length.out = 100)
plot(xs, qgnorm(xs, 0, sqrt(2), 2), type = "l", xlab = "p", ylab = expression(paste("q: p = Pr(", x<=q, ")", sep = "")))

Exponential Power Random Variate Generate, rgnorm

The function rgnorm generates exponential power random variates using the same scale-sign stochastic representation of an exponential power random variable used to compute the inverse CDF. Recall that an exponential power distributed random variable, \(x\), with parameters \(\mu\), \(\alpha\) and \(\beta\) can be written as \(x \stackrel{d}{=} su + \mu\), where \(s\) and \(u\) are independent random variables with \(p(s | \alpha, \beta) = \frac{\beta}{\alpha \Gamma(1/\beta)}\text{exp}\{-(\frac{s}{\alpha})^\beta\}\) and \(u\) uniformly distributed on \(\{-1, 1\}\)

Drawing a value of \(u\) that is uniformly distributed on \(\{-1, 1\}\) is simple. We can draw a value \(s\) according to \(p(s | \alpha, \beta)\) using the inverse CDF function of \(s\). For fixed \(p\), the inverse CDF function of \(s\) finds the value \(q\) that satisfies \(p = \text{Pr}(s \leq q | \alpha, \beta)\):

\[ \begin{aligned} p &= \int_0^q \frac{\beta}{\alpha \Gamma(1/\beta)}\text{exp}\{-(\frac{s}{\alpha})^\beta\}ds \\ &= \int_0^{q^\beta} \frac{1}{\alpha \Gamma(1/\beta)}t^{1/\beta - 1}\text{exp}\{-(\frac{1}{\alpha})^\beta t\}dt & s^\beta = t. \end{aligned} \] As we’ve seen a few times before, this is a gamma CDF. Let \(g(p; \frac{1}{\beta}, (\frac{1}{\alpha})^\beta)\) refer to the inverse CDF function for a gamma distribution with shape \(\frac{1}{\beta}\) and rate \((\frac{1}{\alpha})^\beta\) evaluated at \(p\). Then we can generate \(s\) as follows by drawing a value of \(p\) from a uniform distribution on \((0, 1)\) and setting \(s = g(p; \frac{1}{\beta}, (\frac{1}{\alpha})^\beta)\). This is how rgnorm generates draws from the exponential power distribution. The function rgnorm is a parametrized like rnorm. The first argument is an integer n, which gives the number of draws to take, and the second through fourth arguments are values of \(\mu\), \(\alpha\) and \(\beta\).

Below, we give an example of using the rgnorm to draw exponential power random variates for \(\mu = 0\), \(\alpha = \sqrt{2}\) and \(\beta = 2\). As in the previous examples - this is the standard normal CDF.

xs <- rgnorm(100, 0, sqrt(2), 2)
hist(xs, xlab = "x", freq = FALSE, main = "Histogram of Draws")