1En|y for One-Stage Dorfman-Sterrett

 /****************************************************************/
 /*          S A S   S A M P L E   L I B R A R Y                 */
 /*                                                              */
 /*    NAME: IEDS1EN                                             */
 /*   TITLE: 1En|y for One-Stage Dorfman-Sterrett                */
 /* PRODUCT: QC                                                  */
 /*  SYSTEM: ALL                                                 */
 /*    KEYS: Inspection Sampling,                                */
 /*   PROCS: TABULATE                                            */
 /*    DATA:                                                     */
 /*                                                              */
 /*    MISC:                                                     */
 /*                                                              */
 /*   NOTES: This program tabulates the conditional expected     */
 /*          number of tests 1En|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 ]                       */
 /*                                                              */
 /*                                                              */
 /*     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 n p0 p0prime p pprime result;

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

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

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

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

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

            /*---initialize condexp = 1en|y---*/
            condexp = 0;

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

               if pprime = 0.0 then condexp = 1 + p0prime * n;
               else do;

                  condexp = 1
                            + pprime
                            + ( n * pprime * p0prime )
                            - p0prime;
                  condexp = condexp
                            - ( ( ( 1 - pprime )**( n - 2 ) )*
                            ( 1 + ( 2 * pprime * p0prime ) -
                            ( pprime + p0prime ) ) );
                  condexp = condexp * p0prime / pprime;
                  condexp = condexp + 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 */

               condexp = 1 + p0*(n - outersum);

               end;  /* finish finding condexp when y > 0 */

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