PC(C) for Curtailed Group Testing-Finite

 /****************************************************************/
 /*          S A S   S A M P L E   L I B R A R Y                 */
 /*                                                              */
 /*    NAME: IECURT3                                             */
 /*   TITLE: PC(C) for Curtailed Group Testing-Finite            */
 /* PRODUCT: QC                                                  */
 /*  SYSTEM: ALL                                                 */
 /*    KEYS: Inspection Sampling,                                */
 /*   PROCS: TABULATE                                            */
 /*    DATA:                                                     */
 /*                                                              */
 /*    MISC:                                                     */
 /*                                                              */
 /*   NOTES: This program tabulates the probability of correct   */
 /*          classification of a conforming item 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.3.   */
 /*                                                              */
 /*          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  = 'PC(C)';

   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 - 1, d );
   lognum = lgamma( d + 1 )
            + lgamma( nlot - d )
            + lgamma( nsample )
            + lgamma( ( nlot - nsample ) + 1 );
   lgnlot = lgamma( nlot );

   /*---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 = 1 to nsample;
               do z2 = 0 to nsample - 1;

                  /*---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 - 1;
                  link cond;
                  fact2  = cprob;

                  /*---conditional probability---*/
                  fact3 = ( 1 - pprime ) *
                          min( ( nsample-z1 )/( nsample-z2 ), 1 );
                  fact3 = fact3 +
                          pprime * max((1+z2-z1) / (z2+1),
                                        0 );

                  innersum = innersum + fact1 * fact2 * fact3;

                  end;
               end;

            /*---hypergeometric probability---*/
            logden  = lgamma( y + 1 )
                      + lgamma( ( d - y ) + 1 )
                      + lgamma( nsample - y )
                      + lgamma( ( nlot-d ) - ( nsample-y ) + 1 )
                      + lgnlot;
            hypprob = exp( lognum - logden );

            outersum = outersum + innersum * hypprob;

            end;  /* finish loop over y values */

         result = 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;

   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 ( q>0 ) & ( 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;

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

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