SEVERITY Procedure

Example 28.7 Defining a Finite Mixture Model That Has a Scale Parameter

(View the complete code for this example.)

A finite mixture model is a stochastic model that postulates that the probability distribution of the data generation process is a mixture of a finite number of probability distributions. For example, when an insurance company analyzes loss data from multiple policies that are underwritten in different geographic regions, some regions might behave similarly, but the distribution that governs some regions might be different from the distribution that governs other regions. Further, it might not be known which regions behave similarly. Also, the larger amounts of losses might follow a different stochastic process from the stochastic process that governs the smaller amounts of losses. It helps to model all policies together in order to pool the data together and exploit any commonalities among the regions, and the use of a finite mixture model can help capture the differences in distributions across regions and ranges of loss amounts.

Formally, if f Subscript i and upper F Subscript i denote the PDF and CDF, respectively, of component distribution i and p Subscript i represents the mixing probability that is associated with component i, then the PDF and CDF of the finite mixture of K distribution components are

StartLayout 1st Row 1st Column f left-parenthesis x semicolon normal upper Theta comma bold p right-parenthesis 2nd Column equals sigma-summation Underscript i equals 1 Overscript upper K Endscripts p Subscript i Baseline f Subscript i Baseline left-parenthesis x semicolon normal upper Theta Subscript i Baseline right-parenthesis 2nd Row 1st Column upper F left-parenthesis x semicolon normal upper Theta comma bold p right-parenthesis 2nd Column equals sigma-summation Underscript i equals 1 Overscript upper K Endscripts p Subscript i Baseline upper F Subscript i Baseline left-parenthesis x semicolon normal upper Theta Subscript i Baseline right-parenthesis EndLayout

where normal upper Theta Subscript i denotes the parameters of component distribution i and normal upper Theta denotes the parameters of the mixture distribution, which is a union of all the normal upper Theta Subscript i parameters. bold p denotes the set of mixing probabilities. All mixing probabilities must add up to 1 (sigma-summation Underscript i equals 1 Overscript upper K Endscripts p Subscript i Baseline equals 1).

You can define the finite mixture of a specific number of components and specific distributions for each of the components by defining the FCMP functions for the PDF and CDF. However, in general, it is not possible to fit a scale regression model by using any finite mixture distribution unless you take special care to ensure that the mixture distribution has a scale parameter. This example provides a formulation of a two-component finite mixture model that has a scale parameter.

To start with, each component distribution must have either a scale parameter or a log-transformed scale parameter. Let theta 1 and theta 2 denote the scale parameters of the first and second components, respectively. Let p 1 equals p be the mixing probability, which makes p 2 equals 1 minus p by using the constraint on bold p. The PDF of the mixture of these two distributions can be written as

f left-parenthesis x semicolon theta 1 comma theta 2 comma normal upper Phi comma p right-parenthesis equals StartFraction p Over theta 1 EndFraction f 1 left-parenthesis StartFraction x Over theta 1 EndFraction semicolon normal upper Phi 1 right-parenthesis plus StartFraction 1 minus p Over theta 2 EndFraction f 2 left-parenthesis StartFraction x Over theta 2 EndFraction semicolon normal upper Phi 2 right-parenthesis

where normal upper Phi 1 and normal upper Phi 2 denote the sets of nonscale parameters of the first and second components, respectively, and normal upper Phi denotes a union of normal upper Phi 1 and normal upper Phi 2. For the mixture to have the scale parameter theta, the PDF must be of the form

f left-parenthesis x semicolon theta comma normal upper Phi prime comma p right-parenthesis equals StartFraction 1 Over theta EndFraction left-parenthesis p f 1 left-parenthesis StartFraction x Over theta EndFraction semicolon normal upper Phi prime 1 right-parenthesis plus left-parenthesis 1 minus p right-parenthesis f 2 left-parenthesis StartFraction x Over theta EndFraction semicolon normal upper Phi prime 2 right-parenthesis right-parenthesis

where normal upper Phi prime, normal upper Phi prime 1, and normal upper Phi prime 2 denote the modified sets of nonscale parameters. One simple way to achieve this is to make theta 1 equals theta 2 equals theta and normal upper Phi prime equals normal upper Phi; that is, you simply equate the scale parameters of both components and keep the set of nonscale parameters unchanged. However, forcing the scale parameters to be equal in both components is restrictive, because the mixture cannot model potential differences in the scales of the two components. A better approach is to tie the scale parameters of the two components by a ratio such that theta 1 equals theta and theta 2 equals rho theta. If the ratio parameter rho is estimated along with the other parameters, then the mixture distribution becomes flexible enough to model the variations across the scale parameters of individual components.

To summarize, the PDF and CDF are of the following form for the two-component mixture that has a scale parameter:

StartLayout 1st Row 1st Column f left-parenthesis x semicolon theta comma rho comma normal upper Phi comma p right-parenthesis 2nd Column equals StartFraction 1 Over theta EndFraction left-parenthesis p f 1 left-parenthesis StartFraction x Over theta EndFraction semicolon normal upper Phi 1 right-parenthesis plus left-parenthesis 1 minus p right-parenthesis f 2 left-parenthesis StartFraction x Over theta EndFraction semicolon rho comma normal upper Phi 2 right-parenthesis right-parenthesis 2nd Row 1st Column upper F left-parenthesis x semicolon theta comma rho comma normal upper Phi comma p right-parenthesis 2nd Column equals p upper F 1 left-parenthesis StartFraction x Over theta EndFraction semicolon normal upper Phi 1 right-parenthesis plus left-parenthesis 1 minus p right-parenthesis upper F 2 left-parenthesis StartFraction x Over theta EndFraction semicolon rho comma normal upper Phi 2 right-parenthesis EndLayout

This can be generalized to a mixture of K components by introducing the upper K minus 1 ratio parameters rho Subscript i that relate the scale parameters of each of the K components to the scale parameter theta of the mixture distribution as follows:

StartLayout 1st Row 1st Column theta 1 2nd Column equals theta 2nd Row 1st Column theta Subscript i 2nd Column equals rho Subscript i Baseline theta semicolon i element-of left-bracket 2 comma upper K right-bracket EndLayout

In order to illustrate this approach, define a mixture of two lognormal distributions by using the following PDF function:

StartLayout 1st Row 1st Column f left-parenthesis x semicolon mu comma sigma 1 comma p 2 comma rho 2 comma sigma 2 right-parenthesis equals 2nd Column StartFraction left-parenthesis 1 minus p 2 right-parenthesis Over sigma 1 x StartRoot 2 pi EndRoot EndFraction exp left-parenthesis StartFraction minus left-parenthesis log left-parenthesis x right-parenthesis minus mu right-parenthesis squared Over 2 sigma 1 squared EndFraction right-parenthesis plus 2nd Row 1st Column Blank 2nd Column StartFraction p 2 Over sigma 2 x StartRoot 2 pi EndRoot EndFraction exp left-parenthesis StartFraction minus left-parenthesis log left-parenthesis x right-parenthesis minus mu minus log left-parenthesis rho 2 right-parenthesis right-parenthesis squared Over 2 sigma 2 squared EndFraction right-parenthesis EndLayout

You can verify that mu serves as the log of the scale parameter theta (mu equals log left-parenthesis theta right-parenthesis). The following PROC FCMP steps encode this formulation in a distribution named SLOGNMIX2 for use with PROC SEVERITY:

/*- 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 SEVERITY 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;

As shown in previous examples, an important aspect of defining a distribution for use with PROC SEVERITY is the definition of the PARMINIT subroutine that initializes the parameters. For mixture distributions, in general, the parameter initialization is a nontrivial task. For a two-component mixture, some simplifying assumptions make the problem easier to handle. For the initialization of SLOGNMIX2, the initial values of p 2 and rho 2 are fixed at 0.5, and the following two simplifying assumptions are made:

  • The median of the mixture is the average of the medians of the two components:

    upper F Superscript negative 1 Baseline left-parenthesis 0.5 right-parenthesis equals left-parenthesis exp left-parenthesis mu 1 right-parenthesis plus exp left-parenthesis mu 2 right-parenthesis right-parenthesis slash 2 equals exp left-parenthesis mu right-parenthesis left-parenthesis 1 plus rho 2 right-parenthesis slash 2

    Solution of this equation yields the value of mu in terms of rho 2 and the sample median.

  • Each component has the same mean, which implies the following:

    exp left-parenthesis mu plus sigma 1 squared slash 2 right-parenthesis equals exp left-parenthesis mu plus log left-parenthesis rho 2 right-parenthesis plus sigma 2 squared slash 2 right-parenthesis

    If upper X Subscript i represents the random variable of component distribution i and X represents the random variable of the mixture distribution, then the following equation holds for the raw moment of any order k:

    upper E left-bracket upper X Superscript k Baseline right-bracket equals sigma-summation Underscript i equals 1 Overscript upper K Endscripts p Subscript i Baseline upper E left-bracket upper X Subscript i Superscript k Baseline right-bracket

    This, in conjunction with the assumption on component means, leads to the equations

    StartLayout 1st Row 1st Column log left-parenthesis m 1 right-parenthesis 2nd Column equals mu plus StartFraction sigma 1 squared Over 2 EndFraction 2nd Row 1st Column log left-parenthesis m 1 right-parenthesis 2nd Column equals mu plus log left-parenthesis rho 2 right-parenthesis plus StartFraction sigma 2 squared Over 2 EndFraction EndLayout

    where m 1 denotes the first raw moment of the sample. Solving these equations leads to the following values of sigma 1 and sigma 2:

    StartLayout 1st Row 1st Column sigma 1 squared 2nd Column equals 2 left-parenthesis log left-parenthesis m 1 right-parenthesis minus mu right-parenthesis 2nd Row 1st Column sigma 2 squared 2nd Column equals 2 left-parenthesis log left-parenthesis m 1 right-parenthesis minus mu minus log left-parenthesis rho 2 right-parenthesis right-parenthesis EndLayout

    Note that sigma 1 has a valid value only if log left-parenthesis m 1 right-parenthesis greater-than mu. Among the many possible methods of ensuring this condition, the SLOGNMIX2_PARMINIT subroutine uses the method of doing a linear search over rho 2.

Even when the preceding assumptions are not true for a given problem, they produce reasonable initial values to help guide the nonlinear optimizer to an acceptable optimum if the mixture of two lognormal distributions is indeed a good fit for your input data. This is illustrated by the results of the following steps that fit the SLOGNMIX2 distribution to simulated data, which have different means for the two components (12.18 and 22.76, respectively), and the median of the sample (15.94) is not equal to the average of the medians of the two components (7.39 and 20.09, respectively):

/*-------- 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 severity data=testlognmix print=all plots(histogram kernel)=all;
   loss y;
   dist slognmix2 logn;
run;

The comparison of the fit statistics of SLOGNMIX2 and LOGN, as shown in Output 28.7.1, confirms that the two-component mixture is certainly a better fit to these data than the single lognormal distribution.

Output 28.7.1: Comparison of Fitting One versus Two Lognormal Components to Mixture Data

All Fit Statistics
Distribution-2 Log
Likelihood
AICAICCBICKSADCvM
slognmix238343*38353*38353*38386*0.52221*0.19843*0.02728*
Logn39073 39077 39077 39090 5.86522 66.93414 11.72703 
Note: The asterisk (*) marks the best model according to each column's criterion.


The comparative plot of probability densities in Output 28.7.2 shows that the density function of the mixture distribution is bimodal. In fact, one of the key motivations for using mixture distributions is to find better-fitting models for multimodal data.

The P-P and Q-Q plots in Figure 28 visually confirm that SLOGNMIX2 fits the data very well.

Output 28.7.2: Comparison of PDF Estimates of the Fitted Models

Comparison of PDF Estimates of the Fitted Models


Figure 28: P-P and Q-Q Plots to Evaluate SLOGNMIX2 Fit

sevex07o2g sevex07o2g1


The detailed results for the SLOGNMIX2 distribution are shown in Output 28.7.3. According to the "Initial Parameter Values and Bounds" table, the initial value of rho 2 is not 0.5, indicating that a linear search was conducted to ensure log left-parenthesis m 1 right-parenthesis greater-than mu.

Output 28.7.3: Detailed Estimation Results for the SLOGNMIX2 Distribution

The SEVERITY Procedure
slognmix2 Distribution

Distribution Information
Nameslognmix2
DescriptionMixture of two lognormals with a log-scale parameter Mu
Distribution Parameters5

Initial Parameter Values and Bounds
ParameterInitial
Value
Lower
Bound
Upper
Bound
Mu2.92006-InftyInfty
Sigma10.104551.05367E-8Infty
P20.500001.05367E-81.00000
Rho20.720001.05367E-81.00000
Sigma20.817281.05367E-8Infty

Convergence Status
Convergence criterion (GCONV=1E-8) satisfied.

Optimization Summary
Optimization TechniqueTrust Region
Iterations7
Function Calls18
Log Likelihood-19171.5

Parameter Estimates
ParameterDFEstimateStandard
Error
t ValueApprox
Pr > |t|
Mu13.009220.01554193.68<.0001
Sigma110.495160.0145134.13<.0001
P210.406190.0260015.62<.0001
Rho210.372120.0203818.26<.0001
Sigma211.000190.0212447.09<.0001


By using the relationship that mu 2 equals mu plus log left-parenthesis rho 2 right-parenthesis, you can see that the final parameter estimates are indeed close to the true parameter values that were used to simulate the input sample.

Last updated: June 19, 2025