MATLAB: CDF of exponential power distribution

cdfexponential power distributiongeneralized error distributiongeneralized normal distribution

Hi, I'm trying to compute the CDF of the exponential power distribution, also known as the generalized error distribution or the symmetric generalized normal distribution. I've taken the formula from http://en.wikipedia.org/wiki/Generalized_normal_distribution (see "Version 1"). I've got the PDF working fine, but I've clearly got the CDF wrong as the output doesn't have 0 and 1 as limits (instead, the limits depend on the value of beta). I've copied my code below, can anyone see where I've gone wrong? Thanks so much, Paul.
function p = gnormcdf(x,mu,sigma,beta) % beta is the 'shape' parameter
arg1 = (abs(x-mu)./sigma).^beta;
arg2 = 1./beta;
num = gammainc(arg1,arg2,'lower');
den = 2.*gamma(1./beta);
p = 0.5+sign(x-mu).*(num./den);

Best Answer

Matlab's gammainc normalizes the incomplete gamma function by dividing by gamma(1/beta).  See its documentation.  The Wikipedia formula you are using assumes the incomplete gamma function is not normalized in this manner and divides by gamma(1/beta) in the formula.  Thus, in your code you have effectively divided by gamma(1/beta) twice.  That would certainly throw off your computations.
Roger Stafford
Related Question