PDF of Estimator of Cpk

 /****************************************************************/
 /*       S A S   S A M P L E   L I B R A R Y                    */
 /*                                                              */
 /*    NAME: CPKDEN                                              */
 /*   TITLE: PDF of Estimator of Cpk                             */
 /* PRODUCT: QC                                                  */
 /*  SYSTEM: ALL                                                 */
 /*    KEYS: Capability Analysis, Capability Indices,            */
 /*   PROCS: IML GPLOT                                           */
 /*    DATA:                                                     */
 /*                                                              */
 /*     REF: Chou, Y. and Owen, D. B. (1989). "On the            */
 /*          Distributions of the Estimated Process Capability   */
 /*          Indices".  Communications in Statistics--Theory and */
 /*          Methods, 18. 4549-4560.                             */
 /*                                                              */
 /*          Guirguis, G. H. and Rodriguez, R. N. (1992).        */
 /*          "Computation of Owen's Q Function With Applications */
 /*          to Process Capability Analysis". Journal of Quality */
 /*          Technology 24. To appear.                           */
 /*                                                              */
 /*   NOTES: This example provides a SAS/IML program to compute  */
 /*          the probability density function of the usual       */
 /*          estimator of the process capability index Cpk.  The */
 /*          expression for the density function was derived by  */
 /*          Chou and Owen (1989), and the computational method  */
 /*          is given by Guirguis and Rodriguez (1992).          */
 /*                                                              */
 /*          The parameters of the density are: ds = d/sigma,    */
 /*          ms = |mu-muo|/sigma and n = sample size.  Here      */
 /*                                                              */
 /*             d     = ( USL - LSL ) / 2                        */
 /*             mu0   = ( USL + LSL ) / 2                        */
 /*                                                              */
 /*             LSL   = lower specification limit                */
 /*             USL   = upper specification limit                */
 /*                                                              */
 /*             mu    = process mean                             */
 /*             sigma = process standard deviation               */
 /*                                                              */
 /*          The program computes the density for the estimator  */
 /*          in standardized form (scaled so that its expected   */
 /*          value is zero and its standard deviation is one.)   */
 /*          This enables one to compare the density with that   */
 /*          of the standard normal distribution.                */
 /*                                                              */
 /*                                                              */
 /*    MISC: Depending on the system that 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;

  /* Assign the parameters and set up title for plot */
  n   = 30;
  cpl = 0.1;
  cpu = 0.1;
  ran = 4.0;
  m   = 200;
  title1 'Probability Density of Cpk Estimator for n=30';

  /***************************************************************/
  /*                                                             */
  /*                       IML MODULES                           */
  /*                                                             */
  /***************************************************************/

  /* Support of Non-Central-t distribution */
  start supnct(df,nc,inf,sup);
     ex  = sqrt(df/2)*nc*exp(lgamma((df-1)/2)-lgamma(df/2));
     va  = ((df/(df-2))*(1+(nc*nc)))-(ex*ex);
     sd  = sqrt(va);
     sup = ex+(4*sd);
     inf = ex-(4*sd);
  finish;

  /* Non-central t distribution */
  start cdfnct(x,df,nc,inf,sup);
     if (x >= sup)      then y = 1;
     else if (x <= inf) then y = 0;
          else               y = probt(x,df,nc);
     return(y);
  finish;

  /* Integrand for positive section of density Cpk estimator */
  start fun(x) global(ny,nn,nu,nl);
     y = exp(((nn-2)*log(x))-(0.5*x*x));
     y = Y*(((x*x)-(nn-1))/(nn-1));
     y = y*(probnorm((ny*x)-nl)-probnorm(-(ny*x)+nu));
     return(y);
  finish;

  /* Density Cpk estimator */
  start dencpk(y,n,cpu,cpl)
        global(ny,nn,nu,nl,ng,r,infl,supl,infu,supu,mu,si);
     na = 3*sqrt(n);
     h  = 0.000001;
     ny = (na/sqrt(n-1))*y;
     nn = n;
     if (y <= 0) then do;
        zu = (cdfnct((-na*y)+h,n-1,-nu,infu,supu)
                            -cdfnct(-na*y,n-1,-nu,infu,supu))/h;
        zl = (cdfnct((na*y)+h,n-1,nl,infl,supl)
                            -cdfnct(na*y,n-1,nl,infl,supl))/h;
        z  = na*(zu+zl);
     end;
     else do;
        int = 0.0||(r/y);
        call quad(z,"fun",int) scale = si peak = 5.*mu;
        z = (z*ng)/y;
     end;
     return(z);
  finish;

 /****************************************************************/
 /*                                                              */
 /*               MAIN PROGRAM                                   */
 /*                                                              */
 /****************************************************************/

  r   = (cpl+cpu)*sqrt(n-1)*(0.5);
  ng  = (n-1)/exp(lgamma((n-1)/2)+(log(2)*((n-3)/2)));
  nu  = 3*sqrt(n)*cpu;
  nl  = 3*sqrt(n)*cpl;
  rig = ran;
  lef = -ran;
  inc = (rig-lef)/(m-1);

  run supnct(n-1,-nu,infu,supu);
  run supnct(n-1,nl,infl,supl);

  do x = lef to rig by inc;
     y = exp(-0.5*x*x)/sqrt(2*2*arsin(1));
     if(x = lef) then u = 0||x||y;
     else             u = u//(0||x||y);
  end;

  ds = (3/2)*(cpl+cpu);
  ms = (3/2)*(cpl-cpu);
  f1 = (1/3)*sqrt((n-1)/2)*gamma((n-2)/2)*(1/gamma((n-1)/2));
  f2 = sqrt(2/n)*(1/gamma(0.5))*exp(-n*0.5*ms*ms);
  f3 = ms*(1-(2*probnorm(-sqrt(n)*ms)));
  ex = f1*(ds-f2-f3);
  va = ((n-1)/(9*(n-3)))*(ds**2-(2*ds*(f2+f3))+ms**2+(1/n));
  va = va-(ex*ex);
  sd = sqrt(va);
  mu = ex;
  si = sd;

  do x = lef to rig by inc;
     z = (sd*x)+ex;
     y = dencpk(z,n,cpu,cpl)*sd;
     u = u//(1||x||y);
  end;

  name = {aux x y};
  create density from u[colname = name];
  append from u;
  close density;

  data density;
     set density;
     length curve $ 15;
     if (aux = 0) then curve = "N(0,1)";
     else              curve = "Cpk Estimator";
  run;

  * Plot the density function;
  symbol1 i=join value=none l=1 w=2;
  symbol2 i=join value=none l=3 w=2;
  axis1 order = (-5.00 to 5.00 by 1.00) label = none;
  axis2 order = ( 0.00 to 0.50 by 0.05) label = none;
  proc gplot data = density;
     plot y*x = curve /
        haxis  = axis1
        vaxis  = axis2
        vminor = 0
        hminor = 0
        frame
        ;
  run;

  goptions reset=all;