Model Definition for Burr Distribution
/*--------------------------------------------------------------
SAS Sample Library
Name: svrtburr.sas
Description: Example Program from SAS/ETS User's Guide,
The SEVERITY Procedure
Title: Model Definition for Burr 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.
Further, if you run this sample without any modification to the
names of functions and subroutines, then PROC FCMP writes
warnings of the following nature to the SAS log:
WARNING: Function <name> was defined in a previous package.
Function <name> as defined in the current program will
be used as default when the package is not specified.
You can ignore such warnings; they appear because the functions
defined in this sample are already defined in the input library
Sashelp.Svrtdist.
--------------------------------------------------------------*/
proc fcmp library=sashelp.svrtdist outlib=sasuser.svrtdist.mymodels;
function BURR_DESCRIPTION() $32;
length model $32;
model = "Burr Distribution";
return(model);
endsub;
function BURR_PARMCOUNT();
return(3);
endsub;
function BURR_LOGPDF(x, Theta, Alpha, Gamma);
/* Theta : Scale */
/* Alpha : Shape1 */
/* Gamma : Shape2 */
if (x >= 2.220446E-16) then do; /* constant('MACEPS') */
z = x/Theta;
l1 = log(Alpha*Gamma/x);
l2 = log(z);
l3 = log1px(z**Gamma);
logpdf = l1 + Gamma*l2 - (Alpha+1) * l3;
return (logpdf);
end;
else do;
lbig = 174.673089; /* constant('LOGBIG') */
lsmall = -180.218266; /* constant('LOGSMALL') */
logpdf = lbig + Gamma*lsmall;
if (logpdf < lsmall) then
return (lsmall);
else if (logpdf > lbig) then
return (lbig);
else
return (logpdf);
end;
endsub;
function BURR_PDF(x, Theta, Alpha, Gamma);
/* Theta : Scale */
/* Alpha : Shape1 */
/* Gamma : Shape2 */
logpdf = BURR_LOGPDF(x, Theta, Alpha, Gamma);
if (missing(logpdf)) then
return (.);
if (logpdf < 174.673089) then /* constant('LOGBIG') */
return ( exp(logpdf) );
else
return ( 7.237005E75 ); /* constant('BIG') */
endsub;
function BURR_LOGCDF(x, Theta, Alpha, Gamma);
/* Theta : Scale */
/* Alpha : Shape1 */
/* Gamma : Shape2 */
if (x >= 2.220446E-16) then do; /* constant('MACEPS') */
c = (1 + (x/Theta)**Gamma)**(-Alpha);
if (c < 1) then
return (log1px(-c));
else
return (-180.218266); /* constant('LOGSMALL') */
end;
else
return (-180.218266); /* constant('LOGSMALL') */
endsub;
function BURR_CDF(x, Theta, Alpha, Gamma);
/* Theta : Scale */
/* Alpha : Shape1 */
/* Gamma : Shape2 */
if (x >= 2.220446E-16) then /* constant('MACEPS') */
return (1- (1 + (x/Theta)**Gamma)**(-Alpha));
else
return (0);
endsub;
function BURR_LOGSDF(x, Theta, Alpha, Gamma);
/* Theta : Scale */
/* Alpha : Shape1 */
/* Gamma : Shape2 */
return (-Alpha * log1px((x/Theta)**Gamma));
endsub;
function BURR_SDF(x, Theta, Alpha, Gamma);
/* Theta : Scale */
/* Alpha : Shape1 */
/* Gamma : Shape2 */
return ((1 + (x/Theta)**Gamma)**(-Alpha));
endsub;
subroutine BURR_PARMINIT(dim, x[*], nx[*], F[*], Ftype,
Theta, Alpha, Gamma);
outargs Theta, Alpha, Gamma;
array m[3] / nosymbols;
/* Assume Gamma=2 and compute Alpha & Theta using Method of Moments */
call svrtutil_rawmoments(dim, x, nx, 3, m);
if (missing(m[1])) then do;
Theta = .;
Alpha = .;
Gamma = .;
end;
else do;
t1 = 2*m[3] - 3*m[1]*m[2];
if (t1 < 2.220446E-16) then do; /* constant('MACEPS') */
Theta = sqrt(m[2]);
Alpha = 2;
Gamma = 2;
end;
else do;
Theta = sqrt(m[2]*m[3]/t1);
Alpha = 1 + (m[3]/t1);
Gamma = 2;
end;
end;
endsub;
function BURR_QUANTILE(p, Theta, Alpha, Gamma);
if (p >= 1-1.0e-14) then return (.);
return (Theta * ((1-p)**(-1/Alpha) - 1)**(1/Gamma));
endsub;
function BURR_MEAN(x, Theta, Alpha, Gamma);
if not(Alpha * Gamma > 1) then
return (.); /* first moment does not exist */
return (Theta*gamma(1 + 1/Gamma)*gamma(Alpha - 1/Gamma)/gamma(Alpha));
endsub;
quit;