The exponential power distribution, also known as the generalized normal distribution, was first described in Subbotin (1923)^[Subbotin, M. T. "On the Law of Frequency of Error." Mat. Sb. 31.2 (1923): 206-301.] and rediscovered as the generalized normal distribution in Nadarajah (2005)^[Nadarajah, Saralees. "A generalized normal distribution." Journal of Applied Statistics 32.7 (2005): 685-694.]. It generalizes the Laplace, normal and uniform distributions and is pretty easy to work with in many ways, so it can be very useful. Accordingly, I've made a little `R`

package called `gnorm`

that provides density, CDF and quantile functions for the exponential power distribution as well as random variate generation that are analogous to `dnorm`

, `pnorm`

, `qnorm`

and `rnorm`

. The package can be installed either via CRAN using `install.packages(gnorm)`

or via Github using the `devtools`

package and the command `install_github("maryclare/gnorm")`

. We can load after installation as follows:

```
library(gnorm)
```

Since we're going to be generating some random variables, we should also set a seed:

set.seed(100)

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

This is not necessarily the most natural way to parametrize the exponential power density. One can replace $\alpha$ and $\beta$ with the standard deviation of $x$, $\sigma$, and a shape parameter $q$: $$ p(x | \mu, \sigma, q) = \frac{q}{2\tau}\sqrt{\frac{\Gamma(3/q)}{\Gamma(1/q)^3}}\text{exp}{-(\frac{\Gamma(3/q)}{\Gamma(1/q)})^{q/2}(\frac{|x - \mu|}{\sigma})^q} $$

This is similar but not identically to the parametrization preferred by Box and Tiao (1973)^[Box, G. E. P. and G. C. Tiao. "Bayesian inference in Statistical Analysis." Addison-Wesley Pub. Co., Reading, Mass (1973).]. That said, the parametrization using $\alpha$ and $\beta$ is the one featured on Wikipedia so realistically, it's the one everyone will find when they look up this distribution, so we'll use the $\alpha$ and $\beta$ parametrization for the rest of the vignette.

`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 = "")))

`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)^[Gupta, A. K. and T. Vargas. "Elliptically Contoured Models in Statistics." Kluwer Academic Publishers (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 = "")))

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

There at least one other approach to generating exponential power random variables for any value of $q$. A uniform scale mixture representation of an exponential power distributed random variable, $x$, was originally given in Walker and Guttierez-Pena (1999)^[Walker, S. G. and E. Guttierez-Pena "Robustifying Bayesian Procedures". Bayesian Statistics 6, (1999): 685-710.] and was reprinted in Choy and Walker (2002)^[Choy, S. T. B. and S. G. Walker. "The extended exponential power distribution and Bayesian robustness." Statistics \& Probability Letters 65.3 (2003): 227-232.]. It can be used to generate an exponential power random variate, $x$, as follows:

- Draw $\gamma \sim \text{gamma}(\text{shape} = 1 + 1/\beta, \text{rate}=2^{-\beta/2})$;
- Set $\delta = \alpha\gamma^{1/\beta}/\sqrt{2}$;
- Draw $x \sim \text{uniform}(\mu -\delta, \mu + \delta)$.

**Any scripts or data that you put into this service are public.**

Embedding an R snippet on your website

Add the following code to your website.

For more information on customizing the embed code, read Embedding Snippets.