SAS Employee
Posts: 89

# Additional Distributions in PROC SEVERITY

From time to time we are asked to support additional distributions from within the SEVERITY procedure. While we continue to make every effort to expand the library of predefined distributions, did you know that PROC SEVERITY gives users the ability to write, store and utilize their own distributions for use in the procedure through PROC FCMP?

The following code is taken from the SAS/ETS user's guide and shows how to define a mixed-tail distribution.

proc fcmp library=sashelp.svrtdist outlib=work.sevexmpl.models;

function LOGNGPD_DESCRIPTION() \$256;

length desc \$256;

desc1 = "Lognormal Body-GPD Tail Distribution.";

desc2 = " Mu, Sigma, and Xi are free parameters.";

desc3 = " Xr and Pn are constant parameters.";

desc = desc1 || desc2 || desc3;

return(desc);

endsub;

function LOGNGPD_SCALETRANSFORM() \$3;

length xform \$3;

xform = "LOG";

return (xform);

endsub;

subroutine LOGNGPD_CONSTANTPARM(Xr,Pn);

endsub;

function LOGNGPD_PDF(x, Mu,Sigma,Xi,Xr,Pn);

cutoff = exp(Mu) * Xr;

p = CDF('LOGN',cutoff, Mu, Sigma);

if (x < cutoff + constant('MACEPS')) then do;

return ((Pn/p)*PDF('LOGN', x, Mu, Sigma));

end;

else do;

gpd_scale = p*((1-Pn)/Pn)/PDF('LOGN', cutoff, Mu, Sigma);

h = (1+Xi*(x-cutoff)/gpd_scale)**(-1-(1/Xi))/gpd_scale;

return ((1-Pn)*h);

end;

endsub;

function LOGNGPD_CDF(x, Mu,Sigma,Xi,Xr,Pn);

cutoff = exp(Mu) * Xr;

p = CDF('LOGN',cutoff, Mu, Sigma);

if (x < cutoff + constant('MACEPS')) then do;

return ((Pn/p)*CDF('LOGN', x, Mu, Sigma));

end;

else do;

gpd_scale = p*((1-Pn)/Pn)/PDF('LOGN', cutoff, Mu, Sigma);

H = 1 - (1 + Xi*((x-cutoff)/gpd_scale))**(-1/Xi);

return (Pn + (1-Pn)*H);

end;

endsub;

subroutine LOGNGPD_PARMINIT(dim,x

• ,nx
• ,F
• ,Ftype,
•   Mu,Sigma,Xi,Xr,Pn);

outargs Mu,Sigma,Xi,Xr,Pn;

array xe[1] / nosymbols;

array nxe[1] / nosymbols;

eps = constant('MACEPS');

Pn = 0.8; /* Set mixing probability */

_status_ = .;

call streaminit(56789);

Xb = svrtutil_hillcutoff(dim, x, 100, 25, _status_);

if (missing(_status_) or _status_ = 1) then

Xb = svrtutil_percentile(Pn, dim, x, F, Ftype);

/* prepare arrays for excess values */

i = 1;

do while (i <= dim and x < Xb+eps);

i = i + 1;

end;

dime = dim-i+1;

call dynamic_array(xe, dime);

call dynamic_array(nxe, dime);

j = 1;

do while(i <= dim);

xe = x - Xb;

nxe = nx;

i = i + 1;

j = j + 1;

end;

/* Initialize lognormal parameters */

call logn_parminit(dim, x, nx, F, Ftype, Mu, Sigma);

if (not(missing(Mu))) then

Xr = Xb/exp(Mu);

else

Xr = .;

/* Initialize GPD's shape parameter using excess values */

call gpd_parminit(dime, xe, nxe, F, Ftype, theta_gpd, Xi);

endsub;

subroutine LOGNGPD_LOWERBOUNDS(Mu,Sigma,Xi,Xr,Pn);

outargs Mu,Sigma,Xi,Xr,Pn;

Mu = .; /* Mu has no lower bound */

Sigma = 0; /* Sigma > 0 */

Xi = 0; /* Xi > 0 */

endsub;

quit;

/*----- Simulate a sample for the mixed-tail distribution -----*/

data testmixdist(keep=y label='Lognormal Body-GPD Tail Sample');

call streaminit(45678);

label y='Response Variable';

N = 100;

Mu = 1.5;

Sigma = 0.25;

Xi = 1.5;

Pn = 0.8;

/* Generate data comprising the lognormal body */

Nbody = N*Pn;

do i=1 to Nbody;

y = exp(Mu) * rand('LOGNORMAL')**Sigma;

output;

end;

/* Generate data comprising the GPD tail */

cutoff = quantile('LOGNORMAL', Pn, Mu, Sigma);

gpd_scale = (1-Pn) / pdf('LOGNORMAL', cutoff, Mu, Sigma);

do i=Nbody+1 to N;

y = cutoff + ((1-rand('UNIFORM'))**(-Xi) - 1)*gpd_scale/Xi;

output;

end;

run;

/*--- Set the search path for functions defined with PROC FCMP ---*/

options cmplib=(work.sevexmpl);

/*-------- Fit LOGNGPD model with PROC SEVERITY --------*/

proc severity data=testmixdist print=all plots(histogram kernel)=all

outest=parmest;

loss y;

dist logngpd burr logn gpd;

run;

/*-------- Compute tail cutoff and tail distribution's scale --------*/

data xb_thetat(keep=x_b theta_t);

set parmest(where=(_MODEL_='logngpd' and _TYPE_='EST'));

x_b = exp(Mu) * Xr;

theta_t = (CDF('LOGN',x_b,Mu,Sigma)/PDF('LOGN',x_b,Mu,Sigma)) *

((1-Pn)/Pn);

run;

proc print data=xb_thetat noobs;

run;

Discussion stats
• 0 replies
• 391 views
• 0 likes
• 1 in conversation