2En|y for Two-Stage Dorfman-Sterrett

 /****************************************************************/
 /*          S A S   S A M P L E   L I B R A R Y                 */
 /*                                                              */
 /*    NAME: IEDS2EN                                             */
 /*   TITLE: 2En|y for Two-Stage Dorfman-Sterrett                */
 /* PRODUCT: QC                                                  */
 /*  SYSTEM: ALL                                                 */
 /*    KEYS: Inspection Sampling,                                */
 /*   PROCS: TABULATE                                            */
 /*    DATA:                                                     */
 /*                                                              */
 /*    MISC:                                                     */
 /*                                                              */
 /*   NOTES: This program tabulates the conditional expected     */
 /*          number of tests 2En|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 ]                       */
 /*                                                              */
 /*                                                              */
 /*     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.1.                                     */
 /*                                                              */
 /*          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)'
         p0      = 'p0'
         p0prime = 'p0'''
         p       = 'p'
         pprime  = 'p'''
         result  = '2ENn|y'
         k       = 'k (stage)'
         _y_     = 'y' ;

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

   /*---loop over y values---*/
   do _y_ = 0 to nsample by 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;

               sum      = 0.0;
               upper_m_ = nsample - 2;

               do _m_=1 to upper_m_;

                  /*---find condexp1 = 1en-m|0---*/
                  n = nsample - _m_;
                  y = 0;
                  link enyfun;

                  sum = sum
                        + pprime * ( n - condexp1 ) *
                          ( 1 - pprime)**( _m_ - 1 );
                  end;

               condexp2 = 1 + nsample * p0prime - p0prime * sum;

               end;

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

               outsum   = 0.0;
               upper_m_ = nsample - 2;

               do _m_=1 to upper_m_;

                  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 pmtny2 = p(m,t|n,y)---*/
                     m = _m_;
                     t = _t_;
                     n = nsample;
                     y = _y_;
                     link pcfun;
                     link pncfun;
                     pmtny2 = pc + pnc;

                     /*---find condexp1 = 1en-m|0---*/
                     n = nsample - _m_;
                     y = _y_ - _t_;
                     link enyfun;

                     /*---increment inner sum---*/
                     insum = insum
                             + pmtny2 * ( n - condexp1 );

                     end;

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

               condexp2 = 1 + nsample * p0 - p0 * outsum;
               end;  /* finish general _y_ */

            result = condexp2;
            output;

            end; /* finish loop over pprime values */

         end; /* finish loop over p values */

      end; /* finish loop over y values */

   return;  /* finish main program */


   /*------------------------------------------------------------*/
   /*                                                            */
   /* Calculate condexp1 = 1En|y                                 */
   /*                                                            */
   /* Inputs are:  y, n, pprime, p, p0prime, p0                  */
   /*                                                            */
   /*------------------------------------------------------------*/
   enyfun:

   /*---initialize---*/
   condexp1 = 0;

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

      if pprime = 0.0 then condexp1 = 1 + p0prime * n;
      else do;
         condexp1 = 1
                    + pprime
                    + ( n * pprime * p0prime )
                    - p0prime;
         condexp1 = condexp1
                    - ( ( ( 1 - pprime )**( n - 2 ) )*
                    ( 1 + ( 2 * pprime * p0prime ) -
                    ( pprime + p0prime ) ) );
         condexp1 = condexp1 * p0prime / pprime;
         condexp1 = condexp1 + 1;

         end;

      end;

   /*---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);
         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
                       + pmtny * ( ( n - m ) * ( 1 - p0) - 1 );

            end;  /* finish loop for inner sum */


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

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

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

         /*---increment outer sum---*/
         outersum = outersum
                    + innersum
                    + pmyny * ( ( n - m ) * ( 1 - p0prime ) - 1 );

         end;  /* finish loop for outer sum */

      condexp1 = 1 + p0 * ( n - outersum );
      end;

   return;  /* finish computing condexp1 */


   /*------------------------------------------------------------*/
   /*                                                            */
   /* 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.3 / rts=5;
run;