E(N) for Curtailed Group Testing-Finite

 /****************************************************************/
 /*          S A S   S A M P L E   L I B R A R Y                 */
 /*                                                              */
 /*    NAME: IECURT1                                             */
 /*   TITLE: E(N) for Curtailed Group Testing-Finite             */
 /* PRODUCT: QC                                                  */
 /*  SYSTEM: ALL                                                 */
 /*    KEYS: Inspection Sampling,                                */
 /*   PROCS: TABULATE                                            */
 /*    DATA:                                                     */
 /*                                                              */
 /*    MISC:                                                     */
 /*                                                              */
 /*   NOTES: This program tabulates the expected number of tests */
 /*          for curtailed group testing from a finite lot under */
 /*          an imperfect inspection model.                      */
 /*                                                              */
 /*          Notation:                                           */
 /*                                                              */
 /*          nlot    = size of lot                               */
 /*          d       = number of nonconforming items in lot      */
 /*                                                              */
 /*          n       = sample size                               */
 /*                                                              */
 /*          p1      = Pr[ for group test, declare NC | at least */
 /*                        one NC item ]                         */
 /*                                                              */
 /*          p1prime = 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.      */
 /*          (1986), Statistical Effects of Imperfect Inspection */
 /*          Sampling:  III. Screening (Group Testing), Journal  */
 /*          of Quality Technology 20, 98-124.  See Table 9.1.   */
 /*                                                              */
 /*          Johnson, N. L., Kotz, S., and Wu, X. (1991).        */
 /*          Inspection Errors for Attributes in Quality         */
 /*          Control.  London:  Chapman & Hall.  See Chapter 8.  */
 /*                                                              */
 /****************************************************************/

data table;

   keep nlot d nsample p1 p1prime p pprime result;

   label nlot    = 'N (lot)'
         d       = 'D'
         nsample = 'n (sample)'
         p1      = 'p1'
         p1prime = 'p1'''
         p       = 'p'
         pprime  = 'p'''
         result  = 'E(Number Tests)';

   format result 6.4 ;

   /*---set main parameters---*/
   nlot    = 100;
   nsample = 5;
   d       = 5;
   p1      = 0.95;
   p1prime = 0.05;

   /*---local variables---*/
   ylower = max( 0, nsample + d - nlot );
   yupper = min( nsample, d );

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

      /*---loop over pprime values---*/
      do pprime = 0.0, 0.025, 0.05, 0.075, 0.10 ;

         /*---outer sum---*/
         outersum = 0.0 ;
         do y = ylower to yupper by 1;

            /*---inner double sum---*/
            innersum = 0.0;
            do z1 = 0 to nsample;
               do z2 = 0 to nsample;

                  /*---first conditional probability factor---*/
                  q      = p1;
                  qprime = p1prime;
                  zval   = z1;
                  yval   = y;
                  dval   = d;
                  nval   = nsample;
                  link cond;
                  fact1  = cprob;

                  *--second conditional probability factor--;
                  q      = p;
                  qprime = pprime;
                  zval   = z2;
                  yval   = y;
                  dval   = d;
                  nval   = nsample;
                  link cond;
                  fact2  = cprob;

                  /*---conditional expectation factor--*/
                  if z2 > z1 then
                     fact3 = (( nsample + 1 )*z1)/(z2+1);
                  else if z1 = z2 then
                     fact3 = z1*(nsample-z1)*((1/(z1+1)) +
                             (1/(1+nsample-z1)));
                  else
                     fact3=((nsample+1)*(nsample-z1))/
                           (nsample+1-z2);

                  innersum = innersum + fact1 * fact2 * fact3;

                  end;
               end;

            /*---hypergeometric probability---*/
            bign_ = nlot;
            litn_ = nsample;
            d_    = d;
            y_    = y;
            link hypergmt;

            outersum = outersum + innersum * hypprob;

            end;  /* finish loop over y values */

         result = 1 + outersum ;
         output;

         end;  /* finish loop over pprime values */

      end;  /* finish loop over p values */

   return;  /* finish main program */


   /*------------------------------------------------------------*/
   /*                                                            */
   /* This module computes the conditional probability           */
   /*                                                            */
   /*    cprob   = Pr[ Z = zval | Y = yval ]                     */
   /*                                                            */
   /* where                                                      */
   /*                                                            */
   /*    zval    = number of items classified as defective       */
   /*    yval    = number of actually defective items in sample  */
   /*                                                            */
   /*    dval    = number of defectives in the lot               */
   /*    nval    = sample size                                   */
   /*    nlot    = lot size                                      */
   /*    q       = Pr[ correctly classifying a defective item ]  */
   /*    qprime  = Pr[ incorrectly classifying a good item ]     */
   /*                                                            */
   /*                                                            */
   /*------------------------------------------------------------*/
   cond:

   /*---initialize result to zero---*/
   cprob = 0.0;
   fuzz  = 0.001;

   /*---set limits for subscript---*/
   lolim = max( 0, nval + dval - nlot );
   uplim = min( nval, dval );

   if ( q = 0 ) & ( qprime = 0 ) then do;

      if ( yval >= lolim ) & ( yval <= uplim) then do;

         if zval = 0 then cprob = 1;

         end;

      end;

   else
   if ( q = 0 ) & ( abs( qprime - 1 ) < fuzz ) then do;

      if ( yval >= lolim ) & ( yval <= uplim ) then do;

         if zval = nval - yval then cprob = 1;

         end;

      end;

   else
   if ( abs( q - 1 ) < fuzz ) & ( qprime = 0 ) then do;

      if ( lolim <= yval ) & ( yval <= uplim ) then do;

         if zval = yval then cprob = 1;

         end;

      end;

   else
   if ( abs( q - 1 ) < fuzz ) & ( qprime > 0 ) & ( qprime < 1 )
   then do;

      if ( lolim <= yval ) & ( yval <= uplim ) then do;

         if ( yval <= zval ) & ( zval <= nval ) then do;

            n_ = nval - yval;
            p_ = qprime;
            k_ = zval - yval;

            link binomial;

            cprob = binprob;

            end;

         end;

      end;

   else
   if ( q>0 ) & ( q<1 ) & ( abs( qprime - 1 ) < fuzz ) then do;

      if ( lolim <= yval ) & ( yval <= uplim ) then do;

         if ( nval - yval <= zval ) & ( zval <= nval ) then do;

            n_ = yval;
            p_ = q;
            k_ = zval - ( nval - yval );
            link binomial;

            cprob = binprob;
            end;

         end;

      end;

   else
   if ( abs(q-1) < fuzz ) & ( abs(qprime-1) < fuzz ) then do;

      if ( lolim <= yval ) & ( yval <= uplim ) then do;

         if zval =nval then cprob = 1;

         end;

      end;

   else
   if ( 0 < q ) & ( q < 1) & ( qprime > 0 ) & ( qprime < 1 )
   then do;

      if ( lolim <= yval ) & ( yval <= uplim ) then do;

         if ( 0 <= zval ) & ( zval <= nval ) then do;

            xlo = max( 0, yval + zval - nval );
            xup = min( yval, zval );

            /*---convolution of binomial distributions---*/
            do xval = xlo to xup by 1;

               p_ = q;
               k_ = xval;
               n_ = yval;
               link binomial;
               f1 = binprob;

               p_ = qprime;
               k_ = zval - xval;
               n_ = nval - yval;
               link binomial;
               f2 = binprob;

               cprob =cprob + f1 * f2;

               end;

            end;

         end;

      end;

   else
   if ( q = 0 ) & ( 0 < qprime ) & ( qprime < 1 ) then do;

      if ( yval >= lolim ) & ( yval <= uplim ) then do;

         if ( zval >= 0 ) & ( zval <= nval - yval ) then do;

            n_ = nval - yval;
            p_ = qprime;
            k_ = zval;
            link binomial;

            cprob=binprob;

            end;

         end;

      end;


   else
   if ( qprime = 0 ) & ( 0 < q ) & ( q < 1 ) then do;

      if ( lolim <= yval ) &  ( yval <= uplim ) then do;

         if ( zval >= 0 ) & ( zval <= yval ) then do;

            n_ = yval;
            p_ = q;
            k_ = zval;
            link binomial;

            cprob = binprob ;

            end;

         end;

      end;

   return;  /* finish cond */


   /*---Compute Binomial Probability---*/
   binomial:

   binprob=0.0;

   if n_ = 0 then do;

      if k_ = 0 then binprob = 1.0 ;

      end;

   else
   if n_ > 0 then do;

      if ( k_ > 0 ) & ( k_ < n_ ) then
         binprob = probbnml( p_, n_, k_ ) -
                   probbnml( p_, n_, k_-1 );

      else
      if k_ = n_ then do;
         if ( p_> 0.0 ) & ( p_ < 1.0 ) then
            binprob = p_**n_;
         else if p_ = 1.0 then
            binprob = 1.0;
         end;

      else
      if k_ = 0 then do;
         if ( p_ > 0.0 ) & ( p_ < 1.0 ) then
            binprob = (1.0 - p_)**n_;
         else if p_ = 0.0 then
            binprob = 1.0;
         end;

      end;

   /*---finish binomial computation---*/
   return;

   /*---Compute Hypergeometric Probability---*/
   hypergmt:

      hypprob = 0 ;
      minarg  = max( 0, litn_ + d_ - bign_ );
      maxarg  = min( litn_, d_ );

      if y_ = minarg then

         hypprob = probhypr( bign_, d_, litn_, y_ );

      else
      if ( minarg < y_ ) & ( y_ <= maxarg ) then

         hypprob = probhypr( bign_, d_, litn_, y_     ) -
                   probhypr( bign_, d_, litn_, y_ - 1 );

   /*---finish hypergeometric computation---*/
   return;

run;

proc sort data=table;
   by nlot d nsample p1 p1prime;

proc tabulate data=table noseps;
   by nlot d nsample p1 p1prime;
   class p pprime;
   var result;
   table p, pprime*result=' '*sum=' '*f=8.4 / rts=7;
run;