Estimating Parameters Using the Cramer-von Mises Estimator
/*----------------------------------------------------------------------
SAS Sample Library
Name: sevex04.sas
Description: Example Program from SAS/ETS User's Guide,
The SEVERITY Procedure
Title: Estimating Parameters Using the Cramer-von Mises Estimator
Product: SAS/ETS Software
Keys: Severity Distribution Modeling
PROC: SEVERITY
Notes:
------------------------------------------------------------------------*/
ods graphics on;
/*------- Define Lognormal Body-GPD Tail Mixed 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[i] < 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[j] = x[i] - Xb;
nxe[j] = nx[i];
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 HPSEVERITY by using -------
-------- the Cramer-von Mises minimum distance estimator -------*/
proc severity data=testmixdist obj=cvmobj print=all plots=pp;
loss y;
dist logngpd burr logn gpd;
* Cramer-von Mises estimator (minimizes the distance *
* between parametric and nonparametric estimates) *;
cvmobj = _cdf_(y);
cvmobj = (cvmobj -_edf_(y))**2;
run;