/*-------------------------------------------------------------- SAS Sample Library Name: hpseve02.sas Description: Example program from SAS High Performance Analytics User's Guide, The HPSEVERITY Procedure Title: Defining the Normal Distribution with a Scale Parameter Product: SAS/HPA 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;