Defining the Normal Distribution with a Scale Parameter
/*--------------------------------------------------------------
SAS Sample Library
Name: hsevex02.sas
Description: Example program from SAS/ETS High Performance
Procedures User's Guide, The HPSEVERITY Procedure
Title: Defining the Normal Distribution with a Scale Parameter
Product: SAS High Performance Econometrics Software
Keys: Severity Distribution Modeling
PROC: HPSEVERITY
Notes:
--------------------------------------------------------------*/
/*-------- Define Normal Distribution With Scale Parameter ----------*/
proc fcmp library=sashelp.svrtdist outlib=work.sevexmpl.models;
function normal_s_pdf(x, Sigma, Alpha);
/* Sigma : Scale & Standard Deviation */
/* Alpha : Scaled mean */
return ( exp(-(x/Sigma - Alpha)**2/2) /
(Sigma * sqrt(2*constant('PI'))) );
endsub;
function normal_s_cdf(x, Sigma, Alpha);
/* Sigma : Scale & Standard Deviation */
/* Alpha : Scaled mean */
z = x/Sigma - Alpha;
return (0.5 + 0.5*erf(z/sqrt(2)));
endsub;
subroutine normal_s_parminit(dim, x[*], nx[*], F[*], Ftype, Sigma, Alpha);
outargs Sigma, Alpha;
array m[2] / nosymbols;
/* Compute estimates by using method of moments */
call svrtutil_rawmoments(dim, x, nx, 2, m);
Sigma = sqrt(m[2] - m[1]**2);
Alpha = m[1]/Sigma;
endsub;
subroutine normal_s_lowerbounds(Sigma, Alpha);
outargs Sigma, Alpha;
Alpha = .; /* Alpha has no lower bound */
Sigma = 0; /* Sigma > 0 */
endsub;
quit;
/*--- Simulate a Normal sample affected by Regressors ---*/
data testnorm_reg(keep=y x1-x5 Sigma);
array x{*} x1-x5;
array b{6} _TEMPORARY_ (1 0.5 . 0.75 -2 1);
call streaminit(34567);
label y='Normal Response Influenced by Regressors';
do n = 1 to 100;
/* simulate regressors */
do i = 1 to dim(x);
x(i) = rand('UNIFORM');
end;
/* make x2 linearly dependent on x1 */
x(2) = 5 * x(1);
/* compute log of the scale parameter */
logSigma = b(1);
do i = 1 to dim(x);
if (i ne 2) then
logSigma = logSigma + b(i+1) * x(i);
end;
Sigma = exp(logSigma);
y = rand('NORMAL', 25, Sigma);
output;
end;
run;
/*--- Set the search path for functions defined with PROC FCMP ---*/
options cmplib=(work.sevexmpl);
/*-------- Fit models with PROC HPSEVERITY --------*/
proc hpseverity data=testnorm_reg print=all;
loss y;
scalemodel x1-x5;
dist Normal_s burr logn pareto weibull;
run;
/*--- Refit and compare models with higher limit on iterations ---*/
proc hpseverity data=testnorm_reg print=all;
loss y;
scalemodel x1 x3-x5;
dist Normal_s burr weibull;
nloptions absfconv=2.0e-5 maxiter=100 maxfunc=500;
run;