## Scale Regression with a Finite Mixture Model

```/*--------------------------------------------------------------

SAS Sample Library

Name: hsevex08.sas
Description: Example program from SAS/ETS High Performance
Procedures User's Guide, The HPSEVERITY Procedure
Title: Scale Regression with a Finite Mixture Model
Product: SAS High Performance Econometrics Software
Keys: Severity Distribution Modeling
PROC: HPSEVERITY
Notes:

--------------------------------------------------------------*/

/*- Define Mixture of 2 Lognormal Distributions with a Log-Scale Parameter -*/
proc fcmp library=sashelp.svrtdist outlib=work.sevexmpl.models;
function slognmix2_description() \$128;
return ("Mixture of two lognormals with a log-scale parameter Mu");
endsub;

function slognmix2_scaletransform() \$8;
return ("LOG");
endsub;

function slognmix2_pdf(x, Mu, Sigma1, p2, Rho2, Sigma2);
Mu1 = Mu;
Mu2 = Mu + log(Rho2);
pdf1 = logn_pdf(x, Mu1, Sigma1);
pdf2 = logn_pdf(x, Mu2, Sigma2);
return ((1-p2)*pdf1 + p2*pdf2);
endsub;

function slognmix2_cdf(x, Mu, Sigma1, p2, Rho2, Sigma2);
Mu1 = Mu;
Mu2 = Mu + log(Rho2);
cdf1 = logn_cdf(x, Mu1, Sigma1);
cdf2 = logn_cdf(x, Mu2, Sigma2);
return ((1-p2)*cdf1 + p2*cdf2);
endsub;

subroutine slognmix2_parminit(dim, x[*], nx[*], F[*], Ftype,
Mu, Sigma1, p2, Rho2, Sigma2);
outargs Mu, Sigma1, p2, Rho2, Sigma2;
array m[1] / nosymbols;
p2 = 0.5;
Rho2 = 0.5;
median = svrtutil_percentile(0.5, dim, x, F, Ftype);
Mu = log(2*median/1.5);
call svrtutil_rawmoments(dim, x, nx, 1, m);
lm1 = log(m[1]);

/* Search Rho2 that makes log(sample mean) > Mu */
do while (lm1 <= Mu and Rho2 < 1);
Rho2 = Rho2 + 0.01;
Mu = log(2*median/(1+Rho2));
end;
if (Rho2 >= 1) then
/* If Mu cannot be decreased enough to make it less
than log(sample mean), then revert to Rho2=0.5.
That will set Sigma1 and possibly Sigma2 to missing.
PROC HPSEVERITY replaces missing initial values with 0.001. */
Mu = log(2*median/1.5);

Sigma1 = sqrt(2.0*(log(m[1])-Mu));
Sigma2 = sqrt(2.0*(log(m[1])-Mu-log(Rho2)));
endsub;

subroutine slognmix2_lowerbounds(Mu, Sigma1, p2, Rho2, Sigma2);
outargs Mu, Sigma1, p2, Rho2, Sigma2;
Mu = .; /* Mu has no lower bound */
Sigma1 = 0; /* Sigma1 > 0 */
p2 = 0; /* p2 > 0 */
Rho2 = 0; /* Rho2 > 0 */
Sigma2 = 0; /* Sigma2 > 0 */
endsub;

subroutine slognmix2_upperbounds(Mu, Sigma1, p2, Rho2, Sigma2);
outargs Mu, Sigma1, p2, Rho2, Sigma2;
Mu = .; /* Mu has no upper bound */
Sigma1 = .; /* Sigma1 has no upper bound */
p2 = 1; /* p2 < 1 */
Rho2 = 1; /* Rho2 < 1 */
Sigma2 = .; /* Sigma2 has no upper bound */
endsub;
quit;

/*-------- Simulate a lognormal mixture sample ----------*/
data testlognmix(keep=y);
call streaminit(12345);
Mu1 = 2;
Sigma1 = 1;
i = 0;
do j=1 to 2000;
y = exp(Mu1) * rand('LOGNORMAL')**Sigma1;
output;
end;
Mu2 = 3;
Sigma2 = 0.5;
do j=1 to 3000;
y = exp(Mu2) * rand('LOGNORMAL')**Sigma2;
output;
end;
run;

/*-- Fit and compare scale regression models with 2-component  --*/
/*-- lognormal mixture and the standard lognormal distribution --*/
options cmplib=(work.sevexmpl);

proc hpseverity data=testlognmix print=all;
loss y;
dist slognmix2 logn;
run;

```