Moments of Estimator of Cpk

 /****************************************************************/
 /*       S A S   S A M P L E   L I B R A R Y                    */
 /*                                                              */
 /*    NAME: CPKMOM                                              */
 /*   TITLE: Moments of Estimator of Cpk                         */
 /* PRODUCT: QC                                                  */
 /*  SYSTEM: ALL                                                 */
 /*    KEYS: Capability Analysis, Capability Indices,            */
 /*   PROCS: TABULATE                                            */
 /*    DATA:                                                     */
 /*                                                              */
 /*     REF: Bissell, A. F. (1990).  "How Reliable is Your       */
 /*          Capability Index?".  Applied Statistics 30, 331-    */
 /*          340.                                                */
 /*                                                              */
 /*          W. L. Pearn, S. Kotz, N. L. Johnson (1992).         */
 /*          "Distributional and Inferential Properties of       */
 /*          Process Capability Indices".  Journal of Quality    */
 /*          Technology 24, pp. 216-231.                         */
 /*                                                              */
 /*   NOTES: This program computes the expected value and        */
 /*          standard deviation for the estimator of Cpk.        */
 /*                                                              */
 /*    MISC:                                                     */
 /*                                                              */
 /****************************************************************/

   options ps=60 ls=100;

   data cpk;
     keep ds n ms ex va st;
     label n  = 'n = '
           ds = 'd/sigma'
           ms = '|mu-M|/sigma'
           ex = '  Exp  '
           st = '  Std  '
           va = '  Var  ';
     do n = 30 to 60 by 10;
       f = n-1;
       do ms = 0.0 to 2.0 by 0.5;
         do ds = 2 to 6;
           f1 = (1/3)*(gamma((f-1)/2)/gamma(f/2))*sqrt((f/2));
           f2 = sqrt(2/n)*(1/gamma(0.5))*exp(-n*ms*ms*(1/2));
           f3 = ms*(1-2*probnorm(-ms*sqrt(n)));
           ex = f1*(ds-f2-f3);
           va = (f/(9*(f-2)))*(ds**2-(2*ds*(f2+f3))+ms**2+(1/n));
           va = va-ex**2;
           st = sqrt(va);
           output;
         end;
       end;
     end;
   run;

   proc tabulate data=cpk noseps;
      class n ds ms;
      var ex va st;
      table n, ds, ms*(ex st)*sum =' '*f = 7.4 / rts = 9 condense;
   run;