LCL for Cpk (Chou et al.)


 /****************************************************************/
 /*          S A S   S A M P L E   L I B R A R Y                 */
 /*                                                              */
 /*    NAME: CPKCON1                                             */
 /*   TITLE: LCL for Cpk (Chou et al.)                           */
 /* PRODUCT: QC                                                  */
 /*  SYSTEM: ALL                                                 */
 /*    KEYS: Capability Analysis, Capability Indices,            */
 /*   PROCS: IML                                                 */
 /*    DATA:                                                     */
 /*                                                              */
 /*     REF: Youn-Min Chou, D. B. Owen and S. A. Borrego (1990). */
 /*          "Lower Confidence Limits on Process Capability      */
 /*          Indices".  Journal of Quality Technology 22,        */
 /*          pp. 223-229.  Corrigenda, JQT 24, p. 251.           */
 /*                                                              */
 /*          Kushler, R. H. and Hurley, P. (1992).  "Confidence  */
 /*          Bounds for Capability Indices".  Journal of Quality */
 /*          Technology 24, pp. 188-195.                         */
 /*                                                              */
 /*          Guirguis, G. H. and Rodriguez, R. N. (1992).        */
 /*          "Computation of Owen's Q Function Applied to        */
 /*          Process Capability Analysis".  Journal of Quality   */
 /*          Technology 24, pp. 236-246.                         */
 /*                                                              */
 /*   NOTES: This program computes lower confidence limits for   */
 /*          Cpk using the method of COB (1990).  The program    */
 /*          recreates Table 5 of COB.  Note that the results    */
 /*          can be conservative as discussed by Kushler and     */
 /*          Hurley (1992).  Also note that the limits are       */
 /*          computed as a function of Cpk and n.  The program   */
 /*          can be generalized to compute the limits as a       */
 /*          function of CPL, CPU, and n, following the approach */
 /*          of Guirguis and Rodriguez (1992) (see QC Sample     */
 /*          Program CPKCON2).                                   */
 /*                                                              */
 /*    MISC:                                                     */
 /*                                                              */
 /****************************************************************/

 options ps=60 ls=80;

 proc iml;

 /****************************************************************/
 /*                                                              */
 /*                     ASSIGN PARAMETERS                        */
 /*                                                              */
 /****************************************************************/

 n  = do(30,40,10)`;
 k  = do(0.7,1.0,0.1)`;
 ck = j(nrow(k),nrow(n),0.0);
 p  = 0.95;

 /****************************************************************/
 /*                                                              */
 /*                       IML MODULES                            */
 /*                                                              */
 /****************************************************************/

 /* Integrand */
 start fun(x) global(alpha,delta,df,mlog);
   fx    = probnorm(alpha*x-delta) - probnorm(-alpha*x+delta);
   fx    = fx*exp((df-1)*log(x)-(0.5*x*x)-mlog);
   return(fx);
 finish;

 /* Integral */
 start probd(k,n,ck) global(alpha,delta,df,mlog);
   df    = n - 1;
   alpha = (3.*k*sqrt(n))/sqrt(n-1);
   delta = 3.*ck*sqrt(n);
   int   = 6.*sqrt(n-1)*ck/(6.*k) ;
   peak  = int+1.;
   int   = int ||.p;
   mlog  =   lgamma(0.5*(n-1))+(0.5*(n-3)*log(2));
   call quad( z1, "fun", int  ) peak=peak ;
   return(z1);
 finish;

 /* Root search of probd(k,n,ck) - p */
 start qbd(k,n,p);
    lower1  = 1. / ( 9. * n * k * k ) + .5 /  ( n - 1 )   ;
    lower1  = k*(1-probit( p ) * sqrt(lower1));
    p1      = probd(k,n,lower1);
    lower2  = lower1*(1-0.01/sqrt(n));
 loop:
     p2     = probd(k,n,lower2);
     ratio  = log((1/p-1)/(1/p1-1))/log((1/p1-1)/(1/p2-1))*
              log(lower1/lower2);
     lower  = lower1*exp(ratio);
     lower1 = lower2;
     p1     = p2;
     lower2 = lower;
     if abs(p2-p) > 1.e-5 then go to loop;
   return(lower);
 finish;

 /****************************************************************/
 /*                                                              */
 /*                       MAIN PROGRAM                           */
 /*                                                              */
 /****************************************************************/

 do j=1 to nrow(n);
   do i = 1 to nrow(k);
     ck[i,j] = qbd(k[i],n[j],p);
   end;
 end;

 reset linesize=120 name;
 col =  char(n`,5,0);
 row =  char(k`,5,2);
 mattrib ck rowname=(row)
            colname=(col)
            label={' Cpk \ n = '}
            format=5.3;
 print ck;
 quit;
 run;

 title;