gnorm
The exponential power distribution, also known as the generalized normal distribution, was first described in Subbotin (1923)1 and rediscovered as the generalized normal distribution in Nadarajah (2005)2. 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|μ,α,β)=β2αΓ(1/β)exp{−(|x−μ|α)β}, where μ is the mean of x and α>0 and β>0 are scale and shape parameters.
The function dgnorm
gives the exponential power density given above, p(x|μ,α,β), evaluated at specified values of x, μ, α and β. Below, we give an example of evaluating the exponential power density using dgnorm
for μ=0, α=√2 and β=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, logp(x|μ,α,β) is returned.
This is not necessarily the most natural way to parametrize the exponential power density. One can replace α and β with the standard deviation of x, σ, and a shape parameter q: p(x|μ,σ,q)=q2τ√Γ(3/q)Γ(1/q)3exp{−(Γ(3/q)Γ(1/q))q/2(|x−μ|σ)q}
This is similar but not identically to the parametrization preferred by Box and Tiao (1973)3. That said, the parametrization using α and β 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 α and β parametrization for the rest of the vignette.
pgnorm
The exponential power CDF is derived as follows: Pr(x≤q|μ,α,β)=∫q−∞β2αΓ(1/β)exp{−(|x−μ|α)β}dx=∫q−μ−∞β2αΓ(1/β)exp{−(1α)β|z|β}dz={∫q−μ−∞β2αΓ(1/β)exp{−(1α)β(−z)β}dzq≤μ∫0−∞β2αΓ(1/β)exp{−(1α)β(−z)β+∫q−μ0β2αΓ(1/β)exp{−(1α)βzβ}dzq>μ={∫∞−(q−μ)β2αΓ(1/β)exp{−(1α)βzβ}dzq≤μ∫∞0β2αΓ(1/β)exp{−(1α)βzβdz+∫q−μ0β2αΓ(1/β)exp{−(1α)βzβ}dzq>μ={∫∞0β2αΓ(1/β)exp{−(1α)βzβ}dz−∫−(q−μ)0β2αΓ(1/β)exp{−(1α)βzβ}dzq≤μ∫∞0β2αΓ(1/β)exp{−(1α)βzβ}dz+∫q−μ0β2αΓ(1/β)exp{−(1α)βzβ}dzq>μ=∫∞0β2αΓ(1/β)exp{−(1α)βzβ}dz+sign(q−μ)∫|q−μ|0β2αΓ(1/β)exp{−(1α)βzβ}dz=∫∞0w1/β−12αΓ(1/β)exp{−(1α)βw}dw+sign(q−μ)∫|q−μ|β0w1/β−12αΓ(1/β)exp{−(1α)βw}dwzβ=w=12+sign(q−μ)2∫|q−μ|β0w1/β−1αΓ(1/β)exp{−(1α)βw}dw⏟Gamma CDF, Shape=1β, Rate=(1α)β We can evaluate the exponential power CDF using the gamma CDF - the pgamma
function evaluates the CDF of a gamma distribution with parameters 1β and (1α)β at |q−μ|β.
The function pgnorm
returns the value of the CDF for specified values of q, μ, α and β. 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 Pr(x>q|μ,α,β) when set to FALSE
. The default is lower.tail=TRUE
.Below, we give an example of evaluating the exponential power CDF using pgnorm
for μ=0, α=√2 and β=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)4. We can write xd=su+μ, where s and u are independent random variables with p(s|α,β)=βαΓ(1/β)exp{−(sα)β} and u uniformly distributed on {−1,1}.
Given p, we want to find the value of q that satisfies p=Pr(x≤q|μ,α,β). We can rewrite the problem using s and u as finding p and q that satisfy: |p−0.5|=Pr(s≤|q−μ||α,β)Pr(u=sign(q−μ))=12Pr(s≤|q−μ||α,β)=12∫|q−μ|0βαΓ(1/β)exp{−(sα)β}ds=12∫|q−μ|β01αΓ(1/β)t1/β−1exp{−(1α)βt}dt⏟Gamma CDF, Shape=1β, Rate=(1α)β)sβ=t
Let g(2|p−0.5|;1β,(1α)β) refer to the inverse CDF function for a gamma distribution with shape 1β and rate (1α)β evaluated at 2|p−0.5|. Then the inverse CDF of the exponential power distribution is given by: |q−μ|=g(2|p−0.5|;1β,(1α)β)1/β.
Noting that q>μ when p>0.5 and q≤μ when p≤0.5, we get: q=sign(p−0.5)g(2|p−0.5|;1β,(1α)β)1/β+μ.
This is what is implemented by qgnorm
for specified values of q, μ, α and β. 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 pPr(x>q|μ,α,β) 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 μ=0, α=√2 and β=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 μ, α and β can be written as xd=su+μ, where s and u are independent random variables with p(s|α,β)=βαΓ(1/β)exp{−(sα)β} 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|α,β) using the inverse CDF function of s. For fixed p, the inverse CDF function of s finds the value q that satisfies p=Pr(s≤q|α,β):
p=∫q0βαΓ(1/β)exp{−(sα)β}ds=∫qβ01αΓ(1/β)t1/β−1exp{−(1α)βt}dtsβ=t. As we’ve seen a few times before, this is a gamma CDF. Let g(p;1β,(1α)β) refer to the inverse CDF function for a gamma distribution with shape 1β and rate (1α)β 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;1β,(1α)β). 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 μ, α and β.
Below, we give an example of using the rgnorm
to draw exponential power random variates for μ=0, α=√2 and β=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)5 and was reprinted in Choy and Walker (2002)6. It can be used to generate an exponential power random variate, x, as follows:
Subbotin, M. T. “On the Law of Frequency of Error.” Mat. Sb. 31.2 (1923): 206-301.↩
Nadarajah, Saralees. “A generalized normal distribution.” Journal of Applied Statistics 32.7 (2005): 685-694.↩
Box, G. E. P. and G. C. Tiao. “Bayesian inference in Statistical Analysis.” Addison-Wesley Pub. Co., Reading, Mass (1973).↩
Gupta, A. K. and T. Vargas. “Elliptically Contoured Models in Statistics.” Kluwer Academic Publishers (1993).↩
Walker, S. G. and E. Guttierez-Pena “Robustifying Bayesian Procedures”. Bayesian Statistics 6, (1999): 685-710.↩
Choy, S. T. B. and S. G. Walker. “The extended exponential power distribution and Bayesian robustness.” Statistics & Probability Letters 65.3 (2003): 227-232.↩