Moments of Estimator of Cjkp
/****************************************************************/
/* S A S S A M P L E L I B R A R Y */
/* */
/* NAME: CJKPMOM */
/* TITLE: Moments of Estimator of Cjkp */
/* PRODUCT: QC */
/* SYSTEM: ALL */
/* KEYS: Capability Analysis, Capability Indices, */
/* PROCS: IML */
/* 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 expected value and */
/* variance for the estimator of Cjkp (divided by */
/* Cjkp). */
/* */
/* MISC: */
/* */
/****************************************************************/
options ps=60 ls=80;
proc iml;
start ibeta(x,a,b);
scale = exp(lgamma(a)+lgamma(b)-lgamma(a+b));
y = probbeta(x,a,b)*scale;
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 = (r+1) to (n-r-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;
ban2 = 0.0;
do n = 30 to 50 by 10;
ban1 = 0.0;
do td = 0.4 to 0.8 by 0.2;
ex = expzer(n,td,1);
va = expzer(n,td,2)-(ex*ex);
if (ban1=0) then do;
ban1=1;
ss=ex//va;
end;
else
ss=ss||(ex//va);
end;
if (ban2=0) then do;
ban2=1;
tt=ss;
end;
else
tt=tt//ss;
end;
create cjkp from tt;append from tt;
close cjkp;
reset linesize=120 name;
td = do(0.4,0.8,0.2);
col = char(td,5,2);
row = {
' 30 Exp ',
' Var ',
' 40 Exp ',
' Var ',
' 50 Exp ',
' Var '};
mattrib tt rowname=(row)
colname=(col)
label={'|USL-T|/d'}
format=6.3;
print tt;
quit;
run;