Distribution of Defective Items-Finite

 /****************************************************************/
 /*          S A S   S A M P L E   L I B R A R Y                 */
 /*                                                              */
 /*    NAME: IEDISTZ1                                            */
 /*   TITLE: Distribution of Defective Items-Finite              */
 /* PRODUCT: QC                                                  */
 /*  SYSTEM: ALL                                                 */
 /*    KEYS: Inspection Sampling,                                */
 /*   PROCS: TABULATE                                            */
 /*    DATA:                                                     */
 /*                                                              */
 /*    MISC:                                                     */
 /*                                                              */
 /*   NOTES: This program tabulates the distribution of Z, the   */
 /*          apparent number of nonconforming (defective) items  */
 /*          in a sample drawn from a single, finite lot under   */
 /*          an imperfect inspection model.                      */
 /*                                                              */
 /*          Notation:                                           */
 /*                                                              */
 /*          zprob   = Pr[ Z = z ]                               */
 /*                                                              */
 /*          z       = apparent number of nonconforming items    */
 /*          nlot    = size of lot                               */
 /*          d       = number of nonconforming items in lot      */
 /*          nsample = sample size                               */
 /*          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.      */
 /*          (1985), Statistical Effects of Imperfect Inspection */
 /*          Sampling:  I. Some Basic Distributions, Journal of  */
 /*          Quality Technology 17, 1-31.  See Table 1.          */
 /*                                                              */
 /*          Johnson, N. L., Kotz, S., and Wu, X. (1991).        */
 /*          Inspection Errors for Attributes in Quality         */
 /*          Control.  London:  Chapman & Hall.  See Chapter 2.  */
 /*                                                              */
 /****************************************************************/

data table;

   keep nlot d nsample p pprime z zprob;

   label nlot    = 'N (lot)'
         nsample = 'n (sample)'
         d       = 'D'
         p       = 'p'
         pprime  = 'p'''
         z       = 'z'
         zprob   = 'Pr[ Z = z ]' ;

   format zprob 6.4 ;

   /*---set main parameters---*/
   nlot    = 100;
   d       = 5;
   nsample = 5;

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

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

   /*---loop over range of p values---*/
   do p = 0.75 to 1.00 by 0.05 ;

      /*---loop over range of pprime values---*/
      do pprime = 0.0 to 0.10 by 0.025 ;

         /*---determine Pr[ Z = z ]---*/
         do z = 0 to nsample by 1 ;

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

            /*---output Pr[ Z = z ] ---*/
            output;

            end;  /* finish loop over z */

         end;  /* finish loop over pprime */

      end;  /* finish loop over p */

   /*---finish main program---*/
   return;

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

proc tabulate data=table noseps;
   by nlot d nsample p;
   class z pprime;
   var zprob;
   table z, pprime*zprob=' '*sum=' '*f=8.4 / rts=5;
run;