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;