SIMSYSTEM Procedure

(View the complete code for this example.)

In statistical process control, the three-sigma limits for a random distribution are defined as the mean plus three times the standard deviation and the mean minus three times the standard deviation—symbolically, . Three-sigma limits are very commonly used as a rule of thumb for where the vast majority of values from a process should lie; conversely, they are used to say which outlying values might indicate that the process is not in control. Precisely how much of any given distribution is actually covered by the three-sigma limits varies, although by Chebychev’s inequality this actual coverage is bounded below by . This example is patterned after a study of three-sigma coverage by Wheeler (2010).

The following statements use PROC SIMSYSTEM to compute the theoretical three-sigma coverage for a wide variety of distributions. The key element is the CUMPROB= option in the PROC SIMSYSTEM statement, which computes cumulative probabilities.

data mylib.Grid;
   do s2 = 0 to 10 by 0.05;
      do Kurtosis = 1 to 11 by 0.05;
         Skewness = sqrt(s2);
         output;
         end;
      end;
run;

ods select none;
proc simsystem data=mylib.Grid system=johnson cumprob=-3 3 plot=none;
   ods output Parameters=SKProbJ;
run;
ods select all;

data SKProbJ; set SKProbJ;
   Coverage = P2 - P1;
   s2 = Skewness*Skewness;
   format Coverage percent6.;
run;

The first step creates a fine grid of skewness and kurtosis values. The skewness values are created on a squared scale, to match the results of Wheeler (2010). There are over 40,000 pairs of values in all, although about half of them have and thus lie in the "impossible region" for which no distribution is possible. Still, you probably do not want to display all 20,000 of those values, so the subsequent ODS SELECT NONE and ODS SELECT ALL statements suppress this display for the duration of the PROC SIMSYSTEM run, using an ODS OUTPUT statement to save the distribution parameters to a data set named SKProbJ.

This grid of values is provided as input to PROC SIMSYSTEM. No number of simulations or output data table is specified for the procedure, because you are using the procedure to study the underlying distributions themselves. The SYSTEM=JOHNSON option specifies that this study concerns distributions in the Johnson system; very similar results are obtained using SYSTEM=PEARSON. The CUMPROB= option computes the cumulative probabilities for –3 and 3. Because no value for standard deviation is specified, its value is assumed to be 1 for all distributions, and thus the difference between these two cumulative probabilities is the expected coverage rate. A final DATA step computes this coverage value for each distribution and also recovers the squared skewness values for plotting.

The following statements use PROC SGRENDER to create a plot of the three-sigma coverages computed using the previous statements. They use a predefined graphical template that is specifically designed for plotting the contours of a continuous variable against kurtosis and squared skewness.

proc sgrender data=SKProbJ template=acas.simsystem.Graphics.S2KContourMap;
   dynamic _Title = "Three-Sigma Coverage for Johnson Distributions"
           _XVar  = "s2"
           _YVar  = "Kurtosis"
           _ZVar  = "Coverage";
run;

The resulting plot is shown in Output 31.2.1.

Output 31.2.1: Exploring the Coverage of Three-Sigma Limits

Compare this plot to Figure 3 in Wheeler (2010). As noted there, the coverage is quite high for most of the distributions; only for highly skewed distributions very near the boundary of the "impossible region" does the coverage approach the lower bound of 89%. You can use PROC SIMSYSTEM as follows to see what these distributions look like:

data mylib.LowCoverage;
   do s2 = 7 to 10 by 0.5;
      do Kurtosis = 8 to 11 by 0.5;
         Skewness = sqrt(s2);
         if (1 + s2 < Kurtosis <= 1 + s2 + 0.5) then
            output;
         end;
      end;
run;

proc simsystem data=mylib.LowCoverage system=johnson
               plots=mrmap(skewscale=square);
run;

The DATA step creates another grid of skewness and kurtosis values, which are filtered to lie in the small area of low-coverage distributions shown in Output 31.2.1. The plots that this PROC SIMSYSTEM code produces are shown in Output 31.2.2 and Output 31.2.3.

Output 31.2.2: Low-Coverage Johnson Distribution Map

Output 31.2.3: Low-Coverage Johnson Distribution Densities

Output 31.2.2 confirms that these distributions do indeed lie in the low-coverage region of interest. From Output 31.2.3 you can see that they are all quite extreme. Technically they are "U-shaped" distributions in Johnson’s SB family, but in all cases at least one of the two modes is very sharp. As noted by Wheeler (2010), if your process really does have a distribution like this, then it probably indicates an underlying mixture distribution.