## 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;

```