Model Definition for Gamma Distribution
/*--------------------------------------------------------------
SAS Sample Library
Name: svrtgamm.sas
Description: Example Program from SAS/ETS User's Guide,
The SEVERITY Procedure
Title: Model Definition for Gamma Distribution
Product: SAS/ETS Software
Keys: fitting continuous distributions
PROC: SEVERITY
Notes: If you run this sample program without any modification, then
the Sasuser.Svrtdist library contains functions and subroutines
that are identical to those in the Sashelp.Svrtdist library.
--------------------------------------------------------------*/
proc fcmp outlib=sasuser.svrtdist.models;
function GAMMA_DESCRIPTION()$32;
length model $32;
model = "Gamma Distribution";
return(model);
endsub;
function GAMMA_PARMCOUNT();
return(2);
endsub;
function GAMMA_LOGPDF(x, Theta, Alpha);
/* Theta : Scale */
/* Alpha : Shape */
return(logpdf("GAMMA", x, Alpha, Theta));
endsub;
function GAMMA_PDF(x, Theta, Alpha);
/* Theta : Scale */
/* Alpha : Shape */
return(pdf("GAMMA", x, Alpha, Theta));
endsub;
function GAMMA_CDF(x, Theta, Alpha);
/* Theta : Scale */
/* Alpha : Shape */
eps=3.0e-24; /* relative accuracy */
fpmin=1.0e-37; /* smallest floating point value */
y=x/Theta;
if ((y < 0.0) or (Alpha <= 0.0)) then do;
return(.);
end;
gln = lgamma(Alpha);
if (y < (Alpha+1.0)) then do;
/* begin inline of function gser */
if (y <= 0.0) then do;
return(0.0);
end;
else do;
ap=Alpha;
sum=1.0/Alpha;
del=sum;
do i=1 to 1000;
ap=ap+1;
del=del*(y/ap);
sum=sum+del;
if (abs(del) < abs(sum)*eps) then do;
return(sum*exp(-y+Alpha*log(y)-gln));
end;
end;
end;
return(.);
/* end inline of function gser */
end;
else do;
/* begin inline of function gcf */
b=y+1.0-Alpha;
c=1.0/fpmin;
d=1.0/b;
h=d;
do i= 1 to 1000;
an = -i*(i-Alpha);
b=b+2.0;
d=an*d+b;
if (abs(d) < fpmin) then d=fpmin;
c=b+an/c;
if (abs(c) < fpmin) then c=fpmin;
d=1.0/d;
del=d*c;
h=h*del;
if (abs(del-1.0) < eps) then goto MyBreak;
end;
MyBreak:
return(1.0-(exp(-y+Alpha*log(y)-gln)*h));
/* end inline of function gcf */
end;
endsub;
function GAMMA_LOGCDF(x, Theta, Alpha);
/* Theta : Scale */
/* Alpha : Shape */
cdf = GAMMA_CDF(x, Theta, Alpha);
if (missing(cdf)) then
return (.);
if (cdf > 0) then
return (log(cdf));
else
return (-180.218266); /* constant('LOGSMALL') */
endsub;
function GAMMA_LOGSDF(x, Theta, Alpha);
/* Theta : Scale */
/* Alpha : Shape */
cdf = GAMMA_CDF(x, Theta, Alpha);
if (missing(cdf)) then
return (.);
if (cdf < 1) then
return (log1px(-cdf));
else
return (-180.218266); /* constant('LOGSMALL') */
endsub;
function GAMMA_SDF(x, Theta, Alpha);
/* Theta : Scale */
/* Alpha : Shape */
cdf = GAMMA_CDF(x, Theta, Alpha);
if (missing(cdf)) then
return (.);
return (1-cdf);
endsub;
subroutine GAMMA_PARMINIT(dim, x[*], nx[*], F[*], Ftype, Theta, Alpha);
outargs Theta, Alpha;
m1 = 0;
m2 = 0;
lm = 0;
nlm = 0;
n = 0;
do i=1 to dim;
t1 = nx[i] * x[i];
m1 = m1 + t1;
m2 = m2 + (t1 * x[i]);
if (x[i] > 0) then do;
lm = lm + (nx[i] * log(x[i]));
nlm = nlm + nx[i];
end;
n = n + nx[i];
end;
if (n = 0) then do;
Alpha = .;
Theta = .;
end;
else do;
m1 = m1 / n;
usemom = 1;
if (nlm = n) then do;
/* Compute estimates using approximate maximum likelihood */
s = log(m1) - lm/n;
if (s ne 0) then do;
s1 = (s-3)**2 + 24*s;
if (s1 > 0) then do;
Alpha = (3 - s + sqrt(s1))/(12*s);
if (Alpha > 0) then do;
Theta = m1/Alpha;
usemom = 0;
end;
end;
end;
end;
if (usemom = 1) then do;
/* Compute estimates using method of moments */
m2 = m2 / n;
t1 = m2 - m1**2;
if (t1 < 2.220446E-16) then do; /* constant('MACEPS') */
Alpha = 1;
Theta = m1;
end;
else do;
Alpha = m1**2 / t1;
Theta = t1 / m1;
end;
end;
end;
endsub;
function GAMMA_QUANTILE(p, Theta, Alpha);
/* Theta : Scale */
/* Alpha : Shape */
return(quantile("GAMMA", p, Alpha, Theta));
endsub;
function GAMMA_MEAN(x, Theta, Alpha);
return (Theta*Alpha);
endsub;
quit;