2PC(NC)n|y for Two-Stage Dorfman-Sterrett

 /****************************************************************/
 /*          S A S   S A M P L E   L I B R A R Y                 */
 /*                                                              */
 /*    NAME: IEDS2NC                                             */
 /*   TITLE: 2PC(NC)n|y for Two-Stage Dorfman-Sterrett           */
 /* PRODUCT: QC                                                  */
 /*  SYSTEM: ALL                                                 */
 /*    KEYS: Inspection Sampling,                                */
 /*   PROCS: TABULATE                                            */
 /*    DATA:                                                     */
 /*                                                              */
 /*    MISC:                                                     */
 /*                                                              */
 /*   NOTES: This program tabulates the conditional probability  */
 /*          of correct classification 2PC(NC)n|y for the two-   */
 /*          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_ nsample p0 p0prime p pprime result;

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

   /*---set main parameters---*/
   k       = 2;
   nsample = 10;
   p0      = 0.75;
   p0prime = 0.05;

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

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

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

         outsum   = 0.0;
         upper_m_ = nsample - 2;

            do _m_=1 to upper_m_ by 1;

               insum    = 0.0;
               lower_t_ = max(0,_y_+_m_-nsample);
               upper_t_ = min(_m_,_y_);

               do _t_=lower_t_ to upper_t_ by 1;

                  /* find pnc(m,t|n,y) */
                  m = _m_;
                  t = _t_;
                  n = nsample;
                  y = _y_;
                  link pncfun;
                  pnc2=pnc;

                  /*---special case: _t_=_y_---*/
                  if _t_ = _y_ then

                     /* increment inner sum */
                     insum = insum
                           + pnc2;

                  /*---general case---*/
                  else do;

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

                     /*---find prbpcnc1 = 1pc(nc)n-m|y-t---*/
                     n = nsample - _m_;
                     y = _y_ - _t_;
                     link pcnc1fun;

                     /*---increment inner sum---*/
                     insum = insum
                             + pnc2
                             + (pc2 + pnc2)*(_y_-_t_)*prbpcnc1;
                     end;

                  end;

               /* increment outer sum */
               outsum = outsum + insum;
               end;

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

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

            /*---find pc(n-1,y-1|n,y)---*/
            link pcfun;
            pc2 = pc;
            _term3_ = (pnc + pc)*p;

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

            prbpcnc2 = p0 *
               (outsum+_term1_+_term2_+_term3_+_term4_) / _y_;

            result = prbpcnc2;
            output;

            end; /* finish loop over pprime values */

         end; /* finish loop over p values */

      end; /* finish loop over _y_ values */

   return;  /* finish main program */

   /*------------------------------------------------------------*/
   /*                                                            */
   /* Calculate prbpcnc1 = 1PC(NC)n|y for y = 1, ..., n          */
   /*                                                            */
   /*------------------------------------------------------------*/
   pcnc1fun:

   /*---initialize---*/
   prbpcnc1 = 0.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;

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

   return;  /* finish computing prbpcnc1 */

   /*------------------------------------------------------------*/
   /*                                                            */
   /* 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 nsample p0 p0prime p _y_ pprime;

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