Approximate Confidence Limits for Cpk

/****************************************************************/
/*       S A S   S A M P L E   L I B R A R Y                    */
/*                                                              */
/*    NAME: CPKCON3                                             */
/*   TITLE: Approximate Confidence Limits for Cpk               */
/* PRODUCT: QC                                                  */
/*  SYSTEM: ALL                                                 */
/*    KEYS: Capability Analysis, Capability Indices,            */
/*   PROCS: CAPABILITY                                          */
/*    DATA:                                                     */
/*                                                              */
/* SUPPORT: saswgr                                              */
/*     REF: Bissell, A. F. (1990).  "How Reliable is Your       */
/*          Capability Index?".  Applied Statistics 30, 331-340.*/
/*                                                              */
/*          Kushler, R. H. and Hurley, P. (1992).  "Confidence  */
/*          Bounds for Capability Indices".  Journal of Quality */
/*          Technology 24, pp. 188-195.                         */
/*                                                              */
/*          Zhang, N. F., Stenback, G. A., Wardrop, D. M.       */
/*          (1990).  "Interval Estimation of Process Capability */
/*          Index Cpk".  Comm. in Stat. Theory and Methods 19,  */
/*          4455-4470.                                          */
/*                                                              */
/*   NOTES: This program calculates confidence limits for Cpk   */
/*          using three methods:  Bissell (1990), Equation 6 of */
/*          ZSW (1990), and Equation 8 of ZSW (1990).  See      */
/*          Kushler and Hurley (1992) for an assessment of      */
/*          these methods.                                      */
/*                                                              */
/*    MISC:                                                     */
/*                                                              */
/*  Examples in the documentation were created using the        */
/*    statement:  ods listing style=statistical;                */
/*                                                              */
/****************************************************************/

data Titanium;
   label Hardness = 'Hardness Measurement';
   input Hardness @@;
   datalines;
1.38  1.49  1.43  1.60  1.59
1.34  1.44  1.64  1.83  1.57
1.45  1.74  1.61  1.39  1.63
1.73  1.61  1.35  1.51  1.47
1.46  1.41  1.56  1.40  1.58
1.43  1.53  1.53  1.58  1.62
1.58  1.46  1.26  1.57  1.41
1.53  1.36  1.63  1.36  1.66
1.49  1.55  1.67  1.41  1.39
1.75  1.37  1.36  1.86  1.49
;

proc capability data=Titanium noprint;
   var Hardness;
   specs lsl=0.8 usl=2.4;
   output out=Summary
       n    = n
       mean = mean
       std  = std
       lsl  = lsl
       usl  = usl
       cpk  = cpk
       cpklcl = cpklcl
       cpkucl = cpkucl
       cpl  = cpl
       cpu  = cpu ;

data Summary;
   set Summary;
   length Method $ 16;

   Method = "Bissell";
   lcl = cpklcl;
   ucl = cpkucl;
   output;

   * Assign confidence level;
   level = 0.95;
   aux   = probit( 1 - (1-level)/2 );

   Method = "ZSW Equation 6";
   zsw = log(0.5*n-0.5)
         + ( 2*(lgamma(0.5*n-1)-lgamma(0.5*n-0.5)) );
   zsw = sqrt((n-1)/(n-3)-exp(zsw));
   lcl = cpk*(1-aux*zsw);
   ucl = cpk*(1+aux*zsw);
   output;

   Method = "ZSW Equation 8";
   ds = 3*(cpu+cpl)/2;
   ms = 3*(cpl-cpu)/2;
   f1 = (1/3)*sqrt((n-1)/2)*gamma((n-2)/2)*(1/gamma((n-1)/2));
   f2 = sqrt(2/n)*(1/gamma(0.5))*exp(-n*0.5*ms*ms);
   f3 = ms*(1-(2*probnorm(-sqrt(n)*ms)));
   ex = f1*(ds-f2-f3);
   sd = ((n-1)/(9*(n-3)))*(ds**2-(2*ds*(f2+f3))+ms**2+(1/n));
   sd = sd-(ex*ex);
   sd = sqrt(sd);
   lcl = cpk-aux*sd;
   ucl = cpk+aux*sd;
   output;
run;

title "Approximate 95% Confidence Limits for Cpk";
proc print data = Summary noobs;
   var Method lcl cpk ucl;
run;

proc capability data=Titanium specialindices;
   var Hardness;
   specs lsl=0.8 usl=2.4;
run;