Model Definition for Tweedie Distribution with Scale Parameter
/*--------------------------------------------------------------
SAS Sample Library
Name: svrtstwd.sas
Description: Example Program from SAS/ETS User's Guide,
The SEVERITY Procedure
Title: Model Definition for Tweedie Distribution with Scale Parameter
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.models;
function STWEEDIE_DESCRIPTION() $48;
return ("Tweedie Distribution with Scale Parameter");
endsub;
function STWEEDIE_PARMCOUNT();
return (3);
endsub;
function STWEEDIE_LOGPDF(x, theta, lambda, p);
/* translate into equivalent Tweedie parms using the following:
lambda = mu**(2-p)/(phi*(2-p)) --- Poisson mean
alpha = (2-p)/(p-1) --- gamma shape
theta = phi * (p-1) * mu**(p-1)--- gamma scale */
pm1 = p-1;
alpha = (2-p)/pm1;
mu = theta * alpha * lambda;
phi = theta /(pm1 * mu**pm1);
_delta_ = 1.0e-3;
return (logpdf('TWEEDIE', x, p, mu, phi));
endsub;
function STWEEDIE_PDF(x, theta, lambda, p);
/* translate into equivalent Tweedie parms using the following:
lambda = mu**(2-p)/(phi*(2-p)) --- Poisson mean
alpha = (2-p)/(p-1) --- gamma shape
theta = phi * (p-1) * mu**(p-1)--- gamma scale */
pm1 = p-1;
alpha = (2-p)/pm1;
mu = theta * alpha * lambda;
phi = theta /(pm1 * mu**pm1);
_delta_ = 1.0e-3;
return (pdf('TWEEDIE', x, p, mu, phi));
endsub;
function STWEEDIE_LOGCDF(x, theta, lambda, p);
/* translate into equivalent Tweedie parms using the following:
lambda = mu**(2-p)/(phi*(2-p)) --- Poisson mean
alpha = (2-p)/(p-1) --- gamma shape
theta = phi * (p-1) * mu**(p-1)--- gamma scale */
pm1 = p-1;
alpha = (2-p)/pm1;
mu = theta * alpha * lambda;
phi = theta /(pm1 * mu**pm1);
_delta_ = 1.0e-3;
return (logcdf('TWEEDIE', x, p, mu, phi));
endsub;
function STWEEDIE_CDF(x, theta, lambda, p);
/* translate into equivalent Tweedie parms using the following:
lambda = mu**(2-p)/(phi*(2-p)) --- Poisson mean
alpha = (2-p)/(p-1) --- gamma shape
theta = phi * (p-1) * mu**(p-1)--- gamma scale */
pm1 = p-1;
alpha = (2-p)/pm1;
mu = theta * alpha * lambda;
phi = theta /(pm1 * mu**pm1);
_delta_ = 1.0e-3;
return (cdf('TWEEDIE', x, p, mu, phi));
endsub;
function STWEEDIE_LOGSDF(x, theta, lambda, p);
/* translate into equivalent Tweedie parms using the following:
lambda = mu**(2-p)/(phi*(2-p)) --- Poisson mean
alpha = (2-p)/(p-1) --- gamma shape
theta = phi * (p-1) * mu**(p-1)--- gamma scale */
pm1 = p-1;
alpha = (2-p)/pm1;
mu = theta * alpha * lambda;
phi = theta /(pm1 * mu**pm1);
_delta_ = 1.0e-3;
return (logsdf('TWEEDIE', x, p, mu, phi));
endsub;
function STWEEDIE_SDF(x, theta, lambda, p);
/* translate into equivalent Tweedie parms using the following:
lambda = mu**(2-p)/(phi*(2-p)) --- Poisson mean
alpha = (2-p)/(p-1) --- gamma shape
theta = phi * (p-1) * mu**(p-1)--- gamma scale */
pm1 = p-1;
alpha = (2-p)/pm1;
mu = theta * alpha * lambda;
phi = theta /(pm1 * mu**pm1);
_delta_ = 1.0e-3;
return (sdf('TWEEDIE', x, p, mu, phi));
endsub;
subroutine STWEEDIE_LOWERBOUNDS(theta, lambda, p);
outargs theta, lambda, p;
theta = 0;
lambda = 0;
/* restrict p to the range where Tweedie is equivalent to a
compound Poisson */
p = 1;
endsub;
subroutine STWEEDIE_UPPERBOUNDS(theta, lambda, p);
outargs theta, lambda, p;
theta = .;
lambda = .;
/* restrict p to the range where Tweedie is equivalent to a
compound Poisson */
p = 2;
endsub;
subroutine STWEEDIE_PARMINIT(dim, x[*], nx[*], F[*], Ftype,
theta, lambda, p);
outargs theta, lambda, p;
array m[2] / nosymbols;
/* initialize p */
gscale = .;
alpha = .;
call gamma_parminit(dim, x, nx, F, Ftype, gscale, alpha);
if (not(missing(alpha))) then do;
/* use approximation in 1<p<2 range of alpha = (2-p)/(p-1) */
p = (alpha+2)/(alpha+1);
end;
else do;
p = 1.5;
alpha = 1;
end;
/* initialize mu and phi */
call svrtutil_rawmoments(dim, x, nx, 2, m);
mu = m[1]; /* sample mean */
v = m[2] - m[1]*m[1]; /* sample variance */
phi = v/(mu**p); /* use variance = phi*mu**p property of Tweedie */
theta = phi * (p-1) * mu**(p-1);
lambda = mu/(theta*alpha);
endsub;
function STWEEDIE_QUANTILE(prob, theta, lambda, p);
/* translate into equivalent Tweedie parms using the following:
lambda = mu**(2-p)/(phi*(2-p)) --- Poisson mean
alpha = (2-p)/(p-1) --- gamma shape
theta = phi * (p-1) * mu**(p-1)--- gamma scale */
pm1 = p-1;
alpha = (2-p)/pm1;
mu = theta * alpha * lambda;
phi = theta /(pm1 * mu**pm1);
_delta_ = 1.0e-3;
return (quantile('TWEEDIE', prob, p, mu, phi));
endsub;
function STWEEDIE_MEAN(x, Theta, Lambda, P);
return (Theta* Lambda * (2 - P) / (P - 1));
endsub;
quit;