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;