Resources

1PC(NC)n|y for One-Stage Dorfman-Sterrett

 /****************************************************************/
 /*          S A S   S A M P L E   L I B R A R Y                 */
 /*                                                              */
 /*    NAME: IEDS1NC                                             */
 /*   TITLE: 1PC(NC)n|y for One-Stage Dorfman-Sterrett           */
 /* PRODUCT: QC                                                  */
 /*  SYSTEM: ALL                                                 */
 /*    KEYS: Inspection Sampling,                                */
 /*   PROCS: TABULATE                                            */
 /*    DATA:                                                     */
 /*                                                              */
 /*    MISC:                                                     */
 /*                                                              */
 /*   NOTES: This program tabulates the conditional probability  */
 /*          of correct classification 1PC(NC)n|y for the one-   */
 /*          stage Dorfman-Sterrett procedure under an imperfect */
 /*          inspection model.                                   */
 /*                                                              */
 /*          Notation:                                           */
 /*                                                              */
 /*          n       = number of items in group                  */
 /*          y       = number of truly NC items in group         */
 /*                                                              */
 /*          p0      = Pr[ for group test, declare NC | at least */
 /*                        one NC item ]                         */
 /*                                                              */
 /*          p0prime = Pr[ for group test, declare NC | no NC    */
 /*                        items ]                               */
 /*                                                              */
 /*          p       = Pr[ nonconforming item is classified      */
 /*                        as nonconforming ]                    */
 /*          pprime  = Pr[ conforming item is classified as      */
 /*                        nonconforming ]                       */
 /*                                                              */
 /*          Note that p0prime does not enter into the           */
 /*          computations.                                       */
 /*                                                              */
 /*     REF: Johnson, N. L., Kotz, S., and Rodriguez, R. N.      */
 /*          (1990), Statistical Effects of Imperfect Inspection */
 /*          Sampling:  IV. Modified Dorfman Screening Pro-      */
 /*          cedures, Journal of Quality Technology 22, 128-137. */
 /*          See Table 11.2.                                     */
 /*                                                              */
 /*          Johnson, N. L., Kotz, S., and Wu, X. (1991).        */
 /*          Inspection Errors for Attributes in Quality         */
 /*          Control.  London:  Chapman & Hall.  See Chapter 7.  */
 /*                                                              */
 /****************************************************************/

data table;

   keep k y n p0 p0prime p pprime result;

   label n       = 'n (sample)'
         k       = 'k (stage)'
         p0      = 'p0'
         p0prime = 'p0'''
         p       = 'p'
         pprime  = 'p'''
         result  = '1PC(NC)n|y'
         y       = 'y' ;

   /*---set main parameters---*/
   k       = 1;
   n       = 10;
   p0      = 0.75;
   p0prime = 0.05;

   /*---loop over y values---*/
   do y = 1 to n;

      /*---loop over p values---*/
      do p = 0.75, 0.95;

         /*---loop over pprime values---*/
         do pprime = 0.05, 0.25;

            /* find probpcnc = 1pc(nc)n|y */
            probpcnc=0;

            /* set upper limit for m */
            upperm=n-2;

            /* find outer sum */
            outersum=0.0;
            do m=1 to upperm by 1;

               /* find inner sum */
               innersum=0.0;
               lowert=max(0,y+m-n);
               uppert=min(m,y);
               do t=lowert to uppert by 1;

                  /* find pnc = pnc(m,t|n,y) */
                  link pncfun;

                  /* find pc = pc(m,t|n,y) */
                  link pcfun;

                  /* find pmtny = p(m,t|n,y) */
                  pmtny = pnc + pc;

                  /* increment inner sum */
                  innersum = innersum
                           + (pnc + pmtny*(y-t)*p0*p);

                  end;  /* finish loop for inner sum */

               /* increment outer sum */
               outersum = outersum + innersum;

               end;  /* finish loop for outer sum */

            /* find term2 = pnc(n-1,y-1|n,y) */
            m=n-1;
            t=y-1;
            link pncfun;
            term2 = pnc;

            /* find term3 = p(pnc(n-1,y-1|n,y)+pc(n-1,y-1|n,y)) */
            m=n-1;
            t=y-1;
            link pcfun;
            term3 = p*(pnc+pc);

            /* find term4 = pnc(n,y|n,y) */
            m=n;
            t=y;
            link pncfun;
            term4 = pnc;

            /* find term5 = pnc(n-1,y|n,y) */
            m=n-1;
            t=y;
            link pncfun;
            term5 = pnc;

            probpcnc=p0*(outersum+term2+term3+term4+term5)/y;

            result = probpcnc;
            output;

            end; /* finish loop over pprime values */

         end; /* finish loop over p values */

      end; /* finish loop over y values */

   return;  /* finish main program */

   /*------------------------------------------------------------*/
   /*                                                            */
   /* Calculate pnc = PNC( m, t | n, y )                         */
   /*                                                            */
   /*------------------------------------------------------------*/
   pncfun:

   if      t<=0 then pnc = 0.0;
   else if t>m  then pnc = 0.0;
   else do;
      if p = 1.0 then do;
         if t = 1 then do;
            lognum = lgamma( 1 + n - m )
                   + ( m - 1 ) * log( 1 - pprime )
                   + lgamma( y + 1 )
                   + lgamma( 1 + n - y );
            logden = lgamma( y )
                   + lgamma( 1 + ( n - m ) - ( y - 1 ) )
                   + lgamma( 1 + n );
            pnc = exp( lognum - logden );
            end;
         else pnc = 0.0;
         end;
      else do;
         lognum = lgamma( m )
                + lgamma( 1 + n - m )
                + ( t - 1 ) * log( 1 - p )
                + log( p )
                + ( m - t ) * log( 1 - pprime )
                + lgamma( y + 1 )
                + lgamma( 1 + n - y );
         logden = lgamma( t )
                + lgamma( 1 + m - t )
                + lgamma( 1 + y - t )
                + lgamma( 1 + ( n - m ) - ( y - t ) )
                + lgamma( 1 + n );
         pnc = exp( lognum - logden );
         end;
      end;

   return;  /* finish calculating pnc */

   /*------------------------------------------------------------*/
   /*                                                            */
   /* Calculate pc = PC( m, t | n, y )                           */
   /*                                                            */
   /*------------------------------------------------------------*/
   pcfun:

   if      t >= m then pc = 0.0;
   else if t <  0  then pc = 0.0;
   else do;
      if pprime = 0 then pc = 0.0;
      else if p = 1 then do;
         if t = 0 then do;
            lognum = lgamma( 1 + n - m )
                   + ( m - 1 ) * log( 1 - pprime )
                   + log( pprime )
                   + lgamma( y + 1 )
                   + lgamma( 1 + n - y );
            logden = lgamma( 1 + y )
                   + lgamma( 1 + ( n - m ) - y )
                   + lgamma( 1 + n );
            pc = exp( lognum - logden );
            end;
         else pc = 0.0;
         end;
      else do;
         lognum = lgamma( m )
                + lgamma( 1 + n - m )
                + t * log( 1 - p )
                + (m - t - 1 ) * log( 1 - pprime )
                + log( pprime )
                + lgamma( y + 1 )
                + lgamma( 1 + n - y );
         logden = lgamma( 1 + t )
                + lgamma( m - t )
                + lgamma( 1 + y - t )
                + lgamma( 1 + ( n - m ) - ( y - t ) )
                + lgamma( 1 + n );
         pc = exp( lognum - logden );
         end;
      end;
   return;   /* finish calculating pc */

   run;

proc sort data=table;
   by k n p0 p0prime p y pprime;

proc tabulate data=table noseps;
   by k n p0 p0prime p;
   class y pprime;
   var result;
   table y, pprime*result=' '*sum=' '*f=8.4 / rts=5;
run;