LCL for Cpk (Guirguis, Rodriguez)

 /****************************************************************/
 /*          S A S   S A M P L E   L I B R A R Y                 */
 /*                                                              */
 /*    NAME: CPKCON2                                             */
 /*   TITLE: LCL for Cpk (Guirguis, Rodriguez)                   */
 /* PRODUCT: QC                                                  */
 /*  SYSTEM: ALL                                                 */
 /*    KEYS: Capability Analysis, Capability Indices,            */
 /*   PROCS: CAPABILITY 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    */
 /*          is more general than Table 5 of COB in that it      */
 /*          computes the confidence limit as a function of      */
 /*          CPL, CPU, and n, following the approach of Guirguis */
 /*          and Rodriguez (1992).  Note that this method can    */
 /*          yield conservative results as pointed out by        */
 /*          Kushler and Hurley (1992) and Guirguis and          */
 /*          Rodriguez (1992).                                   */
 /*                                                              */
 /*    MISC:                                                     */
 /*                                                              */
 /****************************************************************/

  options ps=60 ls=80;

  data a;
    input x @@;
    cards;
      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=a  noprint;
    spec lsl=0.8 usl=2.4;
    var x;
    output out  = summary
           n    = n
           mean = mean
           std  = std
           cpu  = cpu
           cpl  = cpl
           cpk  = cpk
  ;

  proc iml;

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

 /* Integrand */
 start fun(x) global(a1,a2,b,nn,mlog);
   fx = probnorm(a1*x-b) - probnorm(-a2*x+b);
   fx = fx*exp((nn-2)*log(x)-(0.5*x*x)-mlog);
   return(fx);
 finish;

 /* Integral */
 start probd(cpl,cpu,n,ck) global(a1,a2,b,nn,mlog);
   nn   = n;
   a1   = (3*cpl*sqrt(n))/sqrt(n-1);
   a2   = (3*cpu*sqrt(n))/sqrt(n-1);
   b    = 3*ck*sqrt(n);
   int  = 2*sqrt(n-1)*ck/(cpl+cpu) ;
   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;

 /* Initial value for the root search */
 start ini(cpl,cpu,n,level);
   cko=0.0;
   paux = 1.0;
   aa:
     if(paux < 0.99 | paux < level) then goto bb;
     cko = cko+(3.0/sqrt(n));
     paux=probd(cpl,cpu,n,cko);
     goto aa;
   bb:
     return(cko);
 finish;

 /* Root search */
 start findck(cpl,cpu,n,level);
   ckr = ini(cpl,cpu,n,level);
   count = 1;
   del   = 1;
   do while( abs(del) > 1.e-5 & count < 100 );
     p1 = probd(cpl,cpu,n,ckr);
     p2 = probd(cpl,cpu,n,ckr*1.001);
     if (p1 = . | p2 = .)
       then
         do;
           ckr = .;
           goto out;
         end;
     der = (p2-p1)*1000/ckr;
     del =  p1/der*log(level/p1);
     ckr = ckr + del;
     count = count+1;
   end;
   if ( count >= 100 )
     then
       cr = .;
   out:
     return(ckr);
 finish;

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

  use summary;
  read all var {n,mean,std,cpu,cpl,cpk};
  level = 0.95;

  reset noname;
  ck = findck(cpl,cpu,n,level);

  print ,;
  rowa = {"     Sample Size            ",
          "     Mean                   ",
          "     Standard Deviation     "};
  rowb = {"     Lower Confidence Limit ",
          "     Cpk                    "};

  auxa = n//mean//std;
  auxb = ck//cpk;
  mattrib auxa rowname=(rowa) format=6.2;
  mattrib auxb rowname=(rowb) format=6.2;
  print "   Summary Statistics";
  print auxa;
  print "   95% Lower Confidence Limit for Cpk";
  print auxb;

  quit;
  run;

  title;