Moments of Estimator of Cpm

 /****************************************************************/
 /*       S A S   S A M P L E   L I B R A R Y                    */
 /*                                                              */
 /*    NAME: CPMMOM                                              */
 /*   TITLE: Moments of Estimator of Cpm                         */
 /* PRODUCT: QC                                                  */
 /*  SYSTEM: ALL                                                 */
 /*    KEYS: Capability Analysis, Capability Indices,            */
 /*   PROCS: TABULATE                                            */
 /*    DATA:                                                     */
 /*     REF: 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 Cpm.        */
 /*                                                              */
 /*    MISC:                                                     */
 /*                                                              */
 /****************************************************************/

   options ps=60 ls=100;

   data cpm;
     keep ds n ms ex va st;
     label n  = 'n = '
           f  = 'f'
           ds = 'd/sigma'
           ms = '|mu-M|/sigma'
           ex = '  Exp  '
           st = '  Std  '
           va = '  Var  ';
     do n = 30 to 50 by 10;
       f = n-1;
       do ms = 0.0 to 2.0 by 0.5;
         hlambda = 0.5*n*(ms**2);
         f1 = sqrt((f+1)/2);
         f2 = (f+1);
         do ds = 2 to 6;
           if(hlambda>0.0) then
           do;
             ex = 0.0;
             va = 0.0;
             j = 0;
             cj = 1.;
             ej = 1.;
             do while ( abs(cj) >= 1.e-7*abs(ex) or
                        abs(ej) >= 1.e-7*abs(va) );
               retain ex va;
               aj = (j*log(hlambda))-lgamma(j+1)-hlambda;
               bj = lgamma((f/2)+j)-lgamma(((f+1)/2)+j);
               cj = exp(aj+bj);
               ex = ex+cj;
               dj = 1/(f-1+(2*j));
               ej = exp(aj)*dj;
               va = va+ej;
               j=j+1;
             end;
             ex = f1*ex*(ds/3);
             va = ((ds/3)**2)*f2*va-(ex**2);
             st = sqrt(va);
           output;
           end;
           else
           do;
             ex = f1*(gamma(f/2)/gamma((f+1)/2))*(ds/3);
             va = (((ds/3)**2)*f2*(1/(f-1)))-(ex**2);
             st = sqrt(va);
             output;
           end;
         end;
       end;
     end;
   run;

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