Compute Cpk Using Johnson SB Distribution

 /****************************************************************/
 /*          S A S   S A M P L E   L I B R A R Y                 */
 /*                                                              */
 /*    NAME: JOHNSYS6                                            */
 /*   TITLE: Compute Cpk Using Johnson SB Distribution           */
 /* PRODUCT: QC                                                  */
 /*  SYSTEM: ALL                                                 */
 /*    KEYS: Capability Analysis, Capability Indices,            */
 /*   PROCS: IML CAPABILITY                                      */
 /*    DATA:                                                     */
 /*                                                              */
 /*   NOTES: Computes Cpk based on fitted Johnson curve.         */
 /*    MISC:                                                     */
 /*                                                              */
 /*     REF: Slifker, J. F. and Shapiro, S. S. (1980), "The      */
 /*          Johnson System: Selection and Parameter Estimation",*/
 /*          Technometrics, Vol. 22, No. 2, pg. 239-246.         */
 /*                                                              */
 /*          Bowman, K. O. and Shenton, L. R. (1983), "Johnson's */
 /*          System of Distributions", Encyclopedia of           */
 /*          Statistical Sciences, Vol. 4, pg. 303-314.          */
 /*          John Wiley & Sons, Inc.                             */
 /*                                                              */
 /*          Rodriguez, R. N. (1992), "Recent Developments in    */
 /*          Process Capability Analysis", Journal of Quality    */
 /*          Technology, Vol. 24, No. 4, pg. 176-187.            */
 /*                                                              */
 /****************************************************************/
options  ps=60 ls=80;

title 'Johnson SB System of Distributions';

 /*---------------------------------------------------------*/
 /*                                                         */
 /*    Generate 150 observations from Johnson Sb System     */
 /*                                                         */
 /*---------------------------------------------------------*/

data sbdata;
   lambda  = 1;
   eta     = 1;
   gamma   = 1;
   epsilon = 2;
   do i = 1 to 150;
      z = rannor(993131);
      a = exp( (z-gamma)/eta);
      x = lambda / ( 1 + ( 1/a ) ) + epsilon;
      output;
      end;
   keep x;
run;

 /*--------------------------------------------------*/
 /*                                                  */
 /* Estimate Johnson Sb parameters using PROC IML    */
 /*                                                  */
 /* The parameters are estimated using the method of */
 /* matching percentiles as described in Slifker and */
 /* Shenton (1980).                                  */
 /*                                                  */
 /*--------------------------------------------------*/

proc iml;

sort sbdata out=sorted by x;
use sorted;

read all var {x} into x;

nobs=nrow(x);

 /*---  Choose z-value for percentile fit ---*/

zval = .524;

p3   = probnorm(3*zval);
p1   = probnorm(zval);
pm1  = probnorm(-1*zval);
pm3  = probnorm(-3*zval);

i4   = p3*nobs  + .5;
i3   = p1*nobs  + .5;
i2   = pm1*nobs + .5;
i1   = pm3*nobs + .5;

x3z  = x[int(i4)] + mod(i4,1)*(x[int(i4)+1] - x[int(i4)]);
x1z  = x[int(i3)] + mod(i3,1)*(x[int(i3)+1] - x[int(i3)]);
xm1z = x[int(i2)] + mod(i2,1)*(x[int(i2)+1] - x[int(i2)]);
xm3z = x[int(i1)] + mod(i1,1)*(x[int(i1)+1] - x[int(i1)]);

m    = x3z - x1z;
n    = xm1z - xm3z;
p    = x1z - xm1z;

 /*--- Estimate Johnson Sb Parameter ---*/

temp    = .5*sqrt( (1 + p/m)*(1 + p/n) );
eta     = zval / log( temp + sqrt(temp*temp - 1) );

temp    = (p/n-p/m)*sqrt((1+p/m)*(1+p/n)-4)/( 2*(p*p/(m*n)-1));
gamma   = eta*log(temp + sqrt(temp*temp + 1) );

lambda  = p*sqrt( ( (1 + p/m)*(1 + p/n)-2)**2-4)/(p*p/(m*n)-1);
epsilon = (x1z + xm1z)/2-lambda/2 + p*(p/n-p/m)/(2*(p*p/(m*n)-1));

 /*---------------------------------------------*/
 /*                                             */
 /*    Compute Cpk from fitted Johnson curve    */
 /*                                             */
 /*    For method of computing Cpk based on a   */
 /*    fitted distribution, see Rodriguez(1992) */
 /*                                             */
 /*---------------------------------------------*/

lsl    = 2.1 ;
usl    = 2.75;

alpha  = .00135;

z      = probit(.5);
a      = exp( (z-gamma)/eta);
median = lambda / ( 1 + ( 1/a ) ) + epsilon;

z      = probit(alpha);
a      = exp( (z-gamma)/eta);
lowp   = lambda / ( 1 + ( 1/a ) ) + epsilon;

z      = probit(1-alpha);
a      = exp( (z-gamma)/eta);
highp  = lambda / ( 1 + ( 1/a ) ) + epsilon;

cpk    = min( (median-lsl)/(median-lowp) ,
              (usl-median)/(highp-median) );

create parms var{eta gamma lambda epsilon cpk};
append;
close parms;

quit;

 /*---------------------------------------*/
 /*                                       */
 /* Use PROC CAPABILITY to save width of  */
 /* histogram bars in OUTFIT= dataset.    */
 /*                                       */
 /*---------------------------------------*/

proc capability data=sbdata noprint;
   hist x / outfit = out1 normal(noprint) noplot;
   run;

 /*------------------------------------------------------------*/
 /*                                                            */
 /* Create annotate dataset that will superimpose the johnson  */
 /* curve and a table containing the parameter estimates onto  */
 /* a histogram.                                               */
 /*                                                            */
 /*------------------------------------------------------------*/

data anno;
   merge parms out1;
   length function style color cepsilon clambda cgamma ceta
          ccpk $ 8 text $ 20;
   function = 'point';
   color    = 'yellow';
   size     = 2;
   xsys     = '2';
   ysys     = '2';
   when     = 'a';
   pi       = 3.14159;
   constant = (100*_width_/lambda)*(eta/sqrt(2*pi));
   start    = _midpt1_ - .5*_width_;
   finish   = _midptn_ + .5*_width_;

   /*--- Draw Johnson Curve ---*/

   do x = start to finish by .01;
      if ( x > epsilon ) then do;
         a     = (x-epsilon)/lambda;
         fy    = log(a / (1 - a ));
         y     = constant*
                 (1/(a*(1-a)))*exp(-.5*(gamma + eta*fy)**2);
         end;
      else y = 0;
      output;
      function = 'draw';
      end;

   /*--- Format Parameters and Cpk with 3 Decimal Places ---*/

   cepsilon = put(epsilon,6.3); clambda = put(lambda,6.3);
   cgamma   = put(gamma,  6.3); ceta    = put(eta,   6.3);
   ccpk     = put(cpk,    6.3);

   /*--- Draw Table containing parameter estimates ---*/

   function = 'move';     xsys = '1';
                          ysys = '1';
                          x = 100; y = 100;        output;
                          xsys = 'A';
                          ysys = 'A';
                          x = -3;  y = -3;         output;
   function = 'cntl2txt';                          output;
   function = 'push';                              output;
   function = 'label';    hsys     = 'A';
                          size     =  1;
                          position = 'a';
                          style    = 'none';
                          color    = 'white';
                          x = 0; y = 0;
                          text = 'Johnson (Sb)';   output;

   function = 'move';     x = -13;   y = .5;       output;
   function = 'draw';     color = 'yellow';
                          size  =  2;
                          x = -4;   y = 0;         output;

   function = 'cntl2txt';                          output;

   function = 'label';    size     =  1;
                          color    = 'white';
                          position = 'c';
               x = 0;  y = -3.5;text = 'Epsilon';  output;
               x = -7; y = -2;  text = 'Lambda';   output;
               x = -6; y = -2;  text = 'Gamma';    output;
               x = -5; y = -2;  text = 'Eta';      output;
               x = -3; y = -2;  text = 'Cpk';      output;

   function = 'pop';                               output;
   function = 'label'; position = 'a';  x = 0;
                       y = -3; text = cepsilon;    output;
                       y = -2; text = clambda;     output;
                       y = -2; text = cgamma;      output;
                       y = -2; text = ceta;        output;
                       y = -2; text = ccpk;        output;

   function = 'move';  x = -18; y = 2;             output;
   function = 'poly';     style = 'solid';
                          color = 'blue';
                          x = 0;  y = 0;           output;
   function = 'polycont'; x = 19; y = 0;           output;
                          x = 19; y = -14;         output;
                          x = 0;  y = -14;         output;
                          x = 0;  y = 0;           output;

   function = 'move';  x = 0;  y = -3;             output;
   function = 'draw';  color = 'white';
                       x = 19; y = 0;              output;

   function = 'move';     xsys = '1';
                          ysys = '1';
                          x = 100; y = 100;        output;
                          xsys = 'A';
                          ysys = 'A';
                          x = -3;  y = -3;         output;

   function = 'move';     x = -13;   y = .5;       output;
   function = 'draw';     color = 'yellow';
                          size  =  2;
                          x = -4;   y = 0;         output;

   keep function xsys  ysys     hsys  x     y
        text     size  position style color when;
   run;

 /*--------------------------------------------------------------*/
 /*                                                              */
 /* Display histogram with Johnson curve and parameter estimates */
 /*                                                              */
 /*--------------------------------------------------------------*/

pattern1 c=blue  v=s;
pattern2 c=red   v=s;
pattern3 c=green v=s;

symbol1  c=white;

options nogstyle;
proc capability data=sbdata annotate=anno noprint;
   spec  lsl = 2.1    usl = 2.75
        clsl = white cusl = white
        llsl = 1     lusl = 2
        wlsl = 2     wusl = 2;
   hist x / cframe = gray
            legend = legend1;
   legend1 cframe=gray cborder=black;
run;
options gstyle;

goptions reset=global;