E(N) for Curtailed Group Testing-Infinite

 /****************************************************************/
 /*          S A S   S A M P L E   L I B R A R Y                 */
 /*                                                              */
 /*    NAME: IECURT1B                                            */
 /*   TITLE: E(N) for Curtailed Group Testing-Infinite           */
 /* 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 an infinite lot    */
 /*          under an imperfect inspection model.                */
 /*                                                              */
 /*          Notation:                                           */
 /*                                                              */
 /*          omega   = proportion of defective 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 10.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 omega nsample p1 p1prime p pprime result;

   label omega   = 'omega'
         nsample = 'n (sample)'
         p1      = 'p1'
         p1prime = 'p1'''
         p       = 'p'
         pprime  = 'p'''
         result  = 'E(Number Tests)';

   format result 6.4 ;

   /*---set main parameters---*/
   omega   = 0.05;
   nsample = 5;
   p1      = 0.95;
   p1prime = 0.05;

   /*---local variables---*/
   ylower = 0;
   yupper = nsample;

   /*---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;
                  nval   = nsample;
                  link cond;
                  fact1  = cprob;

                  *--second conditional probability factor--;
                  q      = p;
                  qprime = pprime;
                  zval   = z2;
                  yval   = y;
                  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;

            /*---binomial probability---*/
            p_    = omega;
            n_    = nsample;
            k_    = y;
            link binomial;

            outersum = outersum + innersum * binprob;

            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  */
   /*                                                            */
   /*    nval    = sample 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;

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

      if zval = 0 then cprob = 1;

      end;

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

      if zval = nval - yval then cprob = 1;

      end;

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

      if zval = yval then cprob = 1;

      end;

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

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

      link binomial;

      cprob = binprob;
      end;

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

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

      cprob = binprob;
      end;

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

      if zval =nval then cprob = 1;
      end;

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

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

         /*---convolution of binomial distributions---*/
         do xval = 0 to zval 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;

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

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

      cprob=binprob;
      end;

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

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

      cprob = binprob ;
      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;

run;

proc sort data=table;
   by omega nsample p1 p1prime;

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