Mean Square Error of Estimator for Cpm

 /****************************************************************/
 /*       S A S   S A M P L E   L I B R A R Y                    */
 /*                                                              */
 /*    NAME: CPMMSE                                              */
 /*   TITLE: Mean Square Error of Estimator for Cpm              */
 /* PRODUCT: QC                                                  */
 /*  SYSTEM: ALL                                                 */
 /*    KEYS: Capability Analysis, Capability Indices,            */
 /*   PROCS: G3D                                                 */
 /*    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 calculates the MSE for the estimator   */
 /*          of Cpm using results of Pearn et al. (1992).  The   */
 /*          MSE is computed as a function of the standardized   */
 /*          parameters                                          */
 /*                     d = ( USL - LSL ) / 2 sigma              */
 /*                     | Mu - M | / sigma                       */
 /*          where M=(USL+LSL)/2.                                */
 /*                                                              */
 /*    MISC:                                                     */
 /*                                                              */
 /****************************************************************/

   options ps=60 ls=80;

   data a;
     label ds  = 'd/Sigma'
           ms  = '|Mu-M|/Sigma'
           mse = 'MSE';
       n = 30;
       f = n-1;
       do ms = 0.0 to 2.0 by 0.1;
         lamda = n*(ms**2);
         f1 = sqrt((f+1)/2)*exp(-lamda/2);
         f2 = (f+1)*exp(-lamda/2);
         do ds = 2 to 6 by 0.1;
           cpm = 1/(3*sqrt(ds*ds+ms*ms));
           if(lamda>0.0) then
           do;
             ex = 0.0;
             va = 0.0;
             j = 0;
             do until (exp(aj) <0.0001);
               retain ex va;
               aj = (j*log(lamda/2))-lgamma(j+1);
               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);
             mse = ((ex-cpm)**2) + 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);
             mse = ((ex-cpm)**2) + va;
             output;
           end;
         end;
       end;
   run;

 title 'Mean Square Error of Cpm Estimator for n=30';
 proc g3d data = a;
   plot ds*ms = mse /
     xticknum = 5
     yticknum = 5
     zticknum = 6
     zmin     = 0.0
     zmax     = 5.0
     rotate   = 45
     tilt     = 70
     grid
     ;
 run;

 goptions reset=all;