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

 /****************************************************************/
 /*          S A S   S A M P L E   L I B R A R Y                 */
 /*                                                              */
 /*    NAME: IEDS2C                                              */
 /*   TITLE: 2PC(C)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(C)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.3.                                     */
 /*                                                              */
 /*          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(C)n|y'
         _y_     = 'y' ;

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

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

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

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

            /*---special handling for _y_=0---*/
            if _y_ = 0 then do;

               outsum   = 0.0;
               upper_m_ = nsample-2;
               do _m_=1 to upper_m_;

                  /*---find probpcc1 = 1pc(c)n-m|0---*/
                  n = nsample-_m_;
                  y = 0;
                  link prbccfun;

                  outsum = outsum
                         + (1 + (nsample-_m_)*(1-probpcc1))*pprime*
                           (1-pprime)**(_m_-1);
                  end;

               probpcc2 = outsum + 2*pprime*(1-p)**(nsample-2);
               probpcc2 = 1 - p0prime*probpcc2/nsample;
               end;

            /*---general _y_---*/
            else do;

               /*---set upper limit for m---*/
               upper_m_ = nsample - 2;

               /*---find outer sum---*/
               outsum = 0.0;
               do _m_=1 to upper_m_ by 1;

                  /*---find inner sum---*/
                  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 pc2 = pc(m,t|n,y)---*/
                     m = _m_;
                     t = _t_;
                     n = nsample;
                     y = _y_;
                     link pcfun;
                     pc2 = pc;

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

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

                     /*---find probpcc1 = 1pc(c)n-m|y-t---*/
                     n = nsample-_m_;
                     y = _y_-_t_;
                     link prbccfun;

                     /*---increment inner sum---*/
                     insum = insum
                           + pc2
                           + pmtny2 *
                             ( nsample + _t_ - ( _m_+_y_ ) ) *
                             ( 1 - probpcc1 );

                     end;  /* finish loop for inner sum */

                  /*---increment outer sum---*/
                  outsum = outsum + insum;
                  end;  /* finish loop for outer sum */

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

               /*---find _term3_ = p'p(n-1,y|n,y)---*/
               link pncfun;
               _term3_ = pprime*(pc + pnc);

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

               /*---find _term5_ = pc(n-1,y-1|n,y)---*/
               m = nsample-1;
               t = _y_-1;
               n = nsample;
               y = _y_;
               link pcfun;
               _term5_ = pc;

               probpcc2  = 1 -
                  p0*(outsum+_term2_+_term3_+_term4_+_term5_)/
                               (nsample-_y_);

               end; /* finish for general y */

            result = probpcc2;
            output;

            end; /* finish loop over pprime values */

         end; /* finish loop over p values */

      end; /* finish loop over _y_ values */

   return;  /* finish main program */

   /*------------------------------------------------------------*/
   /*                                                            */
   /* Calculate 1PC(C)n|y for y = 0, ..., n - 1                  */
   /*                                                            */
   /*------------------------------------------------------------*/
   prbccfun:

   /*---initialize---*/
   probpcc1 = 0.0 ;

   /*---special handling if y=0---*/
   if y = 0 then do;

      probpcc1 = 1 + n*pprime*p0prime - p0prime;
      probpcc1 = probpcc1
                 - (1-p0prime)*(1-2*pprime)*(1-pprime)**(n-2);
      probpcc1 = 1 - p0prime*probpcc1/n;

      end;

   else if y = n then probpcc1 = 0.0 ;

   /*---general case: y>0---*/
   else do;

      /*---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-1);
         do t=lowert to uppert by 1;

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

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

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

            /*---increment inner sum---*/
            innersum = innersum
                     + (pc + pmtny*((n+t)-(m+y))*p0*pprime);

            end;  /* finish loop for inner sum */

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

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

         /*---increment outer sum---*/
         outersum = outersum
                  + innersum
                  + pc
                  + (pc + pnc)*(n-m)*p0prime*pprime;

         end;  /* finish loop for outer sum */

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

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

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

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

      probpcc1=1 -
         p0*(outersum + term2 + term3 + term4 + term5)/(n-y);

      end;

   return;  /* finish computing 1PC(C)n|y */

   /*------------------------------------------------------------*/
   /*                                                            */
   /* 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;