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

```