Acceptance Probabilities for Partial Link Sampling

 /****************************************************************/
 /*          S A S   S A M P L E   L I B R A R Y                 */
 /*                                                              */
 /*    NAME: IEPLINK                                             */
 /*   TITLE: Acceptance Probabilities for Partial Link Sampling  */
 /* PRODUCT: QC                                                  */
 /*  SYSTEM: ALL                                                 */
 /*    KEYS: Inspection Sampling,                                */
 /*   PROCS: TABULATE                                            */
 /*    DATA:                                                     */
 /*                                                              */
 /*    MISC:                                                     */
 /*                                                              */
 /*   NOTES: This program tabulates the acceptance probability   */
 /*          for partial link sampling from multiple lots under  */
 /*          an imperfect inspection model (partial link samp-   */
 /*          ling in the sense introduced by Harishchandra and   */
 /*          Srivenkataramana, 1982).                            */
 /*                                                              */
 /*          Notation:                                           */
 /*                                                              */
 /*          nlot    = size of lot                               */
 /*          d1      = number of nonconforming items in lot i-1  */
 /*          d2      = number of nonconforming items in lot i    */
 /*                                                              */
 /*          nsample = sample size                               */
 /*                                                              */
 /*          c1      = acceptance number                         */
 /*          c2      = rejection  number                         */
 /*          c2p     = acceptance number                         */
 /*                                                              */
 /*          z1      = number of apparent defectives in sample   */
 /*                    from (i-1)st lot                          */
 /*          z2      = number of apparent defectives in first    */
 /*                    sample from ith lot                       */
 /*          z2p     = number of apparent defectives in second   */
 /*                    sample from ith lot                       */
 /*                                                              */
 /*          p       = Pr[ nonconforming item is classified      */
 /*                        as nonconforming ]                    */
 /*          pprime  = Pr[ conforming item is classified as      */
 /*                        nonconforming ]                       */
 /*                                                              */
 /*          accprob = Pr[ acceptance ]                          */
 /*                                                              */
 /*                                                              */
 /*          General Procedure:                                  */
 /*                                                              */
 /*          Take a random sample of size n from the ith lot.    */
 /*          Record the apparent number Zi of defective items.   */
 /*                                                              */
 /*          (a) If Zi <= c1 then accept.                        */
 /*                                                              */
 /*          (b) If Zi > c1 then reject.                         */
 /*                                                              */
 /*          (c) If c1 < Zi <= c2 then take a second sample of   */
 /*              size nsample from the ith lot, and record the   */
 /*              apparent number Zi of defective items.  If      */
 /*                                                              */
 /*              Zi-1 + Zi + Zi' <= c2p ,                        */
 /*                                                              */
 /*              then accept.                                    */
 /*                                                              */
 /*          (d) If c1 < Zi <= c2 and                            */
 /*                                                              */
 /*              Zi-1 + Zi + Zi' > c2p ,                         */
 /*                                                              */
 /*              then reject.                                    */
 /*                                                              */
 /*          Note that in the program, Zi-1 is denoted by z1,    */
 /*          Zi by z2, and Zi' by z2p.                           */
 /*                                                              */
 /*                                                              */
 /*     REF: Harishchandra, K. and Srivenkataramana, T. (1982).  */
 /*          Link Sampling for Attributes, Communications in     */
 /*          Statistics--Theory and Methods 11, 1855-1868.       */
 /*                                                              */
 /*          Johnson, N. L., Kotz, S., and Rodriguez, R. N.      */
 /*          (1986), Statistical Effects of Imperfect Inspection */
 /*          Sampling:  II. Double Sampling and Link Sampling,   */
 /*          Journal of Quality Technology 18, 116-138.          */
 /*          See Table 5.                                        */
 /*                                                              */
 /*          Johnson, N. L., Kotz, S., and Wu, X. (1991).        */
 /*          Inspection Errors for Attributes in Quality         */
 /*          Control.  London:  Chapman & Hall.  See Chapter 4.  */
 /*                                                              */
 /****************************************************************/

data table;

   keep nlot d1 d2 nsample c1 c2 c2p p pprime accprob;

   label nlot    = 'N (lot)'
         d1      = 'D(i-1)'
         d2      = 'D(i)'
         nsample = 'n'
         c1      = 'c1'
         c2      = 'c2'
         c2p     = 'c2'''
         p       = 'p'
         pprime  = 'p'''
         accprob = 'Pr[ Accept ]';

   format zprob   6.4
          accprob 6.4 ;

   /*---set main parameters---*/
   nlot    = 100;
   nsample = 20;
   d1      = 15;
   d2      = 10;
   c1      = 1;
   c2      = 5;
   c2p     = 5;

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

      /*---loop over pprime values---*/
      do pprime = 0.0, 0.01, 0.02, 0.05, 0.10;

         /*---compute term1---*/
         link first;

         /*---computer term2---*/
         link second;

         accprob = term1 + term2;

         output;

         end;  /* finish loop over pprime values */

      end;  /* finish loop over p values */

   return;  /* finish main program */

   /*------------------------------------------------------------*/
   /*                                                            */
   /* This module computes the probability Pr[ Z2 <= c1 ]        */
   /*                                                            */
   /* The following serve as input parameters:                   */
   /*                                                            */
   /*    c1      = acceptance value for first sample             */
   /*    nlot    = lot size                                      */
   /*    d       = number of defectives in lot                   */
   /*    nsample = sample size                                   */
   /*    p       = Pr[ correctly classifying a defective item ]  */
   /*    pprime  = Pr[ incorrectly classifying a good item ]     */
   /*                                                            */
   /* The following is returned:                                 */
   /*                                                            */
   /*    term1   = Pr[ Z2 <= c1 ]                                */
   /*                                                            */
   /*------------------------------------------------------------*/
   first:

   term1 = 0.0 ;
   do z2 = 0 to c1 by 1;

      /* nlot, nsample, p, pprime are globally defined */
      z = z2;
      d = d2;
      link uncond;

      term1 = term1 + zprob;

      end;

   return;  /* finish first */

   /*------------------------------------------------------------*/
   /*                                                            */
   /* This module computes the probability                       */
   /*                                                            */
   /*    term2   = Pr[ c1 < Z2 <= c2, Z1 + Z2 + Z2p <= c2p ]     */
   /*                                                            */
   /*                                                            */
   /* The following serve as input parameters:                   */
   /*                                                            */
   /*    c2      = acceptance value                              */
   /*    c2p     = acceptance value                              */
   /*    nlot    = lot size                                      */
   /*    d1      = number of defectives in lot i-1               */
   /*    d2      = number of defectives in lot i                 */
   /*    nsample = sample size                                   */
   /*    p       = Pr[ correctly classifying a defective item ]  */
   /*    pprime  = Pr[ incorrectly classifying a good item ]     */
   /*                                                            */
   /*                                                            */
   /*------------------------------------------------------------*/
   second:

   term2 = 0.0 ;
   zimax = max( c2, c2p );

   do z1 = 0 to zimax by 1;

      do z2 = 0 to zimax by 1;

         do z2p = 0 to zimax by 1;

            if ( c1 < z2 ) & ( z2 <= c2 ) & ( z1+z2+z2p <= c2p )
            then do;

               /* nlot, nsample, p, pprime are globally defined */
               z = z1;
               d = d1;
               link uncond;
               product = zprob;

               /* obtain joint probability that Z2=z2, Z2P=z2p */
               link uncond2;
               product = product * unprb2;

               term2 = term2 + product;
               end;

            end;

         end;

      end;

   return;  /* finish second */

   /*------------------------------------------------------------*/
   /*                                                            */
   /* This module computes the unconditional joint distribution  */
   /* of Z2 and Z2P.                                             */
   /*                                                            */
   /* The following serve as input parameters:                   */
   /*                                                            */
   /*    z2      = number of items classified as defective       */
   /*    z2p     = number of items classified as defective       */
   /*    nlot    = lot size                                      */
   /*    d2      = number of defectives in lot                   */
   /*    nsample = sample size                                   */
   /*    p       = Pr[ correctly classifying a defective item ]  */
   /*    pprime  = Pr[ incorrectly classifying a good item ]     */
   /*                                                            */
   /* The following is returned:                                 */
   /*                                                            */
   /*    unprb2  =                                               */
   /*                                                            */
   /*------------------------------------------------------------*/
   uncond2:

   unprb2 = 0.0 ;
   upp1   = min( d2, nsample );
   lsum   = max( 0, 2 * nsample + d2 - nlot );
   usum   = min( d2, 2 * nsample );

   do ylocal = 0 to upp1 by 1;

      do yplocal = 0 to upp1 by 1;

         if ( lsum <= ylocal + yplocal ) &
            ( ylocal + yplocal <= usum )
         then do;

            /*---absolute hypergeometric probability---*/
            bign_ = nlot;
            litn_ = 2 * nsample;
            d_    = d2;
            y_    = ylocal + yplocal;
            link hypergmt;
            hprob1 = hypprob;

            bign_ = 2 * nsample;
            litn_ = nsample;
            d_    = ylocal + yplocal;
            y_    = ylocal;
            link hypergmt;
            hprob2 = hypprob;

            mhyp = hprob1 * hprob2;

            /*--conditional probability that Z2 = z2 ---*/
            nval = nsample;
            zval = z2;
            yval = ylocal;
            dval = d2;
            link cond;
            mhyp = mhyp * cprob;

            /*---conditional probability that Z2P = z2p ---*/
            nval = nsample;
            zval = z2p;
            yval = yplocal;
            dval = d2;
            link cond;
            mhyp=mhyp * cprob;

            *--add over y and yp--;
            unprb2 = unprb2 + mhyp;

            end;

         end;

      end;

   return;  /* finish uncond2 */


   /*------------------------------------------------------------*/
   /*                                                            */
   /* 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                                      */
   /*    p       = Pr[ correctly classifying a defective item ]  */
   /*    pprime  = Pr[ incorrectly classifying a good item ]     */
   /*                                                            */
   /*                                                            */
   /*------------------------------------------------------------*/
   cond:

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

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

   if ( p = 0 ) & ( pprime = 0 ) then do;

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

         if zval = 0 then cprob = 1;

         end;

      end;

   else
   if ( p = 0 ) & ( abs( pprime - 1 ) < fuzz ) then do;

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

         if zval = nval - yval then cprob = 1;

         end;

      end;

   else
   if ( abs( p - 1 ) < fuzz ) & ( pprime = 0 ) then do;

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

         if zval = yval then cprob = 1;

         end;

      end;

   else
   if ( abs( p - 1 ) < fuzz ) & ( pprime > 0 ) & ( pprime < 1 )
   then do;

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

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

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

            link binomial;

            cprob = binprob;

            end;

         end;

      end;

   else
   if ( p > 0 ) & ( p < 1 ) & ( abs( pprime-1 ) < fuzz ) then do;

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

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

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

            cprob = binprob;
            end;

         end;

      end;

   else
   if ( abs( p-1 ) < fuzz ) & ( abs( pprime-1 ) < fuzz ) then do;

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

         if zval =nval then cprob = 1;

         end;

      end;

   else
   if ( 0 < p ) & ( p < 1) & ( pprime > 0 ) & ( pprime < 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_ = p;
               k_ = xval;
               n_ = yval;
               link binomial;
               f1 = binprob;

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

               cprob =cprob + f1 * f2;

               end;

            end;

         end;

      end;

   else
   if ( p = 0 ) & ( 0 < pprime ) & ( pprime < 1 ) then do;

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

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

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

            cprob=binprob;

            end;

         end;

      end;


   else
   if ( pprime = 0 ) & ( 0 < p ) & ( p < 1 ) then do;

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

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

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

            cprob = binprob ;

            end;

         end;

      end;

   return;  /* finish cond */


   /*------------------------------------------------------------*/
   /*                                                            */
   /* This module computes the probability Pr[ Z = z ], where Z  */
   /* is the number of items classified as defective.            */
   /*                                                            */
   /* The following serve as input parameters:                   */
   /*                                                            */
   /*    z       = number of items classified as defective       */
   /*    nlot    = lot size                                      */
   /*    d       = number of defectives in lot                   */
   /*    nsample = sample size                                   */
   /*    p       = Pr[ correctly classifying a defective item ]  */
   /*    pprime  = Pr[ incorrectly classifying a good item ]     */
   /*                                                            */
   /* The following is returned:                                 */
   /*                                                            */
   /*    zprob   = Pr[ Z = z ]                                   */
   /*                                                            */
   /*------------------------------------------------------------*/
   uncond:

   /*---used for roundoff---*/
   fuzz = 0.0001 ;

   /*---lower and upper limits for y---*/
   miny = max( 0, nsample + d - nlot );
   maxy = min( nsample, d );

   /*---initialize probability to zero---*/
   zprob = 0.0 ;

   /*---Case I: p = 0 ---*/
   if p = 0 then do;

      /*---Ia: pprime = 0 ---*/
      if pprime = 0 then do;

         if z = 0 then zprob = 1 ;

         end;  /* finish Ia */

      /*---Ib: pprime = 1 ---*/
      else if abs( pprime - 1 ) < fuzz then do;

         minz = max( 0, nsample - d );
         maxz = min( nsample, nlot - d );

         if ( minz <= z ) & ( z <= maxz ) then do;

            bign_ = nlot;
            litn_ = nsample;
            d_    = d;
            y_    = nsample - z;
            link hypergmt;

            zprob = hypprob;

            end;

         end;  /* finish Ib */

      /*---Ic:  0 < pprime < 1 ---*/
      else do;

         /* Note: minz =  0 */
         maxz = nsample - max( 0, nsample + d - nlot );

         if ( z <= maxz ) then
         do y = miny to maxy by 1;

            /*---obtain Pr[ Y = y ]---*/
            bign_ = nlot;
            litn_ = nsample;
            d_    = d;
            y_    = y;
            link hypergmt;

            /*---obtain Pr[ Z = z | Y = y ]---*/
            n_ = nsample - y;
            k_ = z;
            p_ = pprime;
            link binomial;

            zprob = zprob + binprob * hypprob ;

            end;

         end;  /* finish Ic */

      end;  /* finish Case I */


   /*---Case II:  p = 1 ---*/
   else if ( abs( p - 1 ) < fuzz ) then do;

      /*---IIa:  pprime = 0 (perfect inspection) ---*/
      if pprime = 0 then do;

         minz = max( 0, nsample + d - nlot );
         maxz = min( nsample, d );

         if ( minz <= z ) & ( z <= maxz ) then do;

            bign_ = nlot;
            litn_ = nsample;
            d_    = d;
            y_    = z;
            link hypergmt;

            zprob = hypprob;

            end;

         end;  /* finish IIa */

      /*---IIb:  pprime = 1 ---*/
      else if ( abs( pprime - 1 ) < fuzz ) then do;

         if z = nsample then zprob = 1 ;

         end;  /* finish IIb */

      /*---IIc:  0 < pprime < 1 ---*/
      else do;

         minz = max( 0, nsample + d - nlot );
         maxz = nsample ;

         if ( minz <= z ) & ( z <= maxz ) then
         do y = miny to maxy by 1;

            /*---compute Pr[ Y = y ] ---*/
            bign_ = nlot ;
            litn_ = nsample ;
            d_    = d;
            y_    = y;
            link hypergmt;

            /*---obtain Pr[ Z = z | Y = y ]---*/
            p_ = pprime;
            k_ = z - y;
            n_ = nsample - y;
            link binomial;

            zprob = zprob + hypprob * binprob;

            end;

         end;  /* finish IIb */

      end;  /* finish Case II */


   /*---Case III:  0 < p < 1---*/
   else do;

      /*---IIIa:  pprime = 0 ---*/
      if pprime = 0 then do;

         /* zmin = 0 */
         zmax = min( nsample, d );

         if z <= zmax then
         do y = miny to maxy by 1;

            /*---obtain Pr[ Y = y ]---*/
            bign_ = nlot ;
            litn_ = nsample ;
            d_    = d;
            y_    = y;
            link hypergmt;

            /*---obtain Pr[ Z = z | Y = y ]---*/
            p_ = p;
            k_ = z;
            n_ = y;
            link binomial;

            /*---increment unconditional probability---*/
            zprob = zprob + binprob * hypprob ;

            end;

         end;  /* finish IIIa */

      /*---IIIb:  pprime = 1 ---*/
      else if abs( pprime - 1 ) < fuzz then do;

         zmin = nsample - min( nsample, d );
         /* zmax = nsample */

         if z >= zmin then
         do y = miny to maxy by 1;

            /*---obtain Pr[ Y = y ]---*/
            bign_ = nlot;
            litn_ = nsample;
            d_    = d;
            y_    = y;
            link hypergmt;

            /*---obtain Pr[ Z = z | Y = y ]---*/
            p_ = p;
            k_ = z - ( nsample - y );
            n_ = y;
            link binomial;

            /*---increment unconditional probability---*/
            zprob = zprob + binprob * hypprob ;

            end;

         end;  /* finish Case IIIb */

      /*---IIIc:  0 < pprime < 1 ---*/
      else
      do y = miny to maxy by 1;

         /*---obtain Pr[ Y = y ]---*/
         bign_ = nlot;
         litn_ = nsample;
         d_    = d;
         y_    = y;
         link hypergmt;

         /*---obtain Pr[ Z = z | Y = y ]---*/
         condprob = 0.0 ;
         minx     = max( 0, y + z - nsample );
         maxx     = min( y, z );

         do x = minx to maxx by 1;

            p_ = p;
            k_ = x;
            n_ = y;
            link binomial;
            factor1 = binprob;

            p_ = pprime;
            k_ = z - x;
            n_ = nsample - y;
            link binomial;
            factor2 = binprob;

            condprob = condprob + factor1 * factor2;

            end;

         /*---increment unconditional probability---*/
         zprob = zprob + condprob * hypprob ;

         end;  /* finish IIIc */

      end;  /* finish Case III */

   return;  /* finish uncond */

   /*---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 d1 d2 nsample c1 c2 c2p;

proc tabulate data=table noseps;
   by nlot d1 d2 nsample c1 c2 c2p;
   class p pprime;
   var accprob;
   table p, pprime*accprob=' '*sum=' '*f=8.4 / rts=7;
run;