Mean Square Error of Cjkp Estimator

 /****************************************************************/
 /*       S A S   S A M P L E   L I B R A R Y                    */
 /*                                                              */
 /*    NAME: CJKPMSE                                             */
 /*   TITLE: Mean Square Error of Cjkp Estimator                 */
 /* PRODUCT: QC                                                  */
 /*  SYSTEM: ALL                                                 */
 /*    KEYS: Capability Analysis, Capability Indices,            */
 /*   PROCS: IML  G3D                                            */
 /*    DATA:                                                     */
 /*                                                              */
 /*     REF: N. L. Johnson, S. Kotz, and W. L. Pearn (1992).     */
 /*          "Flexible Process Capability Indices".  Institute   */
 /*          of Statistics Mimeo Series #2072.  Department of    */
 /*          Statistics, University of North Carolina.           */
 /*                                                              */
 /*   NOTES: This program computes the MSE of the Cjkp estima-   */
 /*          tor.  The MSE is computed as a function of the      */
 /*          parameters ds = d/Sigma and td = (USL-T)/d,  where  */
 /*          d=(USL+LSL)/2.                                      */
 /*                                                              */
 /*    MISC: Depending on the system you use to run this sample, */
 /*          you may find that the run time for this program is  */
 /*          long.                                               */
 /*                                                              */
 /****************************************************************/

options ps=60 ls=80;

proc iml;

  start fun(t) global(aa,bb);
    y = ((aa-1)*log(t))+((bb-1)*log(1-t));
    y = exp(y);
    return(y);
  finish;

  start ibeta(x,a,b) global(aa,bb);
    aa = a;
    bb = b;
    int  = 0.0 || x;
    call quad( y, "fun", int );
    return(y);
  finish;

  start expzer(n,td,r);
    a1 = 1/(td*td);
    y1 = ((0.5*r)*log(n)) +
         lgamma(0.5*(n-r))-((n+r)*log(2))-lgamma(0.5*n);
    y1 = exp(y1);
    y2 = ((2*sqrt(a1))-1 )##r;
    b = 1/((4*a1)-(4*sqrt(a1))+2);
    sum = 0.0;
    do k = 1 to (n-1);
      f1 = lgamma(n+1)+lgamma(0.5*n)-lgamma(n-k+1)-
           lgamma(k+1)-lgamma(0.5*k)-lgamma(0.5*(n-k));
      f1 = exp(f1);
      f2 = ibeta(1-b,0.5*(n-k),0.5*(k-r));
      f3 = ibeta(b,0.5*k,0.5*(n-k-r));
      f4 = f1*(f2+(y2*f3));
      sum = sum+f4;
    end;
    z = y1*(y2+1+sum);
    return(z);
  finish;

  n = 30;
  ban1 = 0.0;
  do td = 0.4 to 1.0 by 0.025;
    e1 = expzer(n,td,1);
    e2 = expzer(n,td,2);
    do ds = 2 to 6 by 0.25;
      cjkp = ds/(3*td);
      mse = (cjkp**2)*( ((e1-1)**2) + (e2-(e1**2)) );
      if (ban1=0) then do;
          ban1=1;
          tt=td||ds||mse;
        end;
      else
        tt=tt//(td||ds||mse);
    end;
  end;
  col = {'td','ds','mse'};
  create a from tt[colname=col];append from tt;
  close a;
quit;

data a;
  set a;
  label td  = '(USL-T)/d'
        ds  = 'd/Sigma'
        mse = 'MSE';

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

goptions reset=all;