Graff-Roeloffs' Modification of Dorfman

 /****************************************************************/
 /*          S A S   S A M P L E   L I B R A R Y                 */
 /*                                                              */
 /*    NAME: IEGRAFF                                             */
 /*   TITLE: Graff-Roeloffs' Modification of Dorfman             */
 /* PRODUCT: QC                                                  */
 /*  SYSTEM: ALL                                                 */
 /*    KEYS: Inspection Sampling,                                */
 /*   PROCS: TABULATE                                            */
 /*    DATA:                                                     */
 /*                                                              */
 /*    MISC:                                                     */
 /*                                                              */
 /*   NOTES: This program tabulates the probabilities of correct */
 /*          classification PC(NC) and PC(C) of nonconforming    */
 /*          and conforming items respectively, and the percent  */
 /*          reduction in expected overall number of tests for   */
 /*          the modified Dorfman procedure proposed by Graff    */
 /*          and Roeloffs.  This consists of repeating inspec-   */
 /*          tion until either r1 NC decisions or r2 C decisions */
 /*          have been obtained, whichever comes first.  Here a  */
 /*          hierarchal version is considered.                   */
 /*                                                              */
 /*          Notation:                                           */
 /*                                                              */
 /*          k       = stage                                     */
 /*                                                              */
 /*          r1j     = r1 for jth stage                          */
 /*          r2j     = r2 for jth stage                          */
 /*                                                              */
 /*          n0      = number of items in group                  */
 /*          n1      = number of items in primary subgroup       */
 /*          n2      = number of items in secondary subgroup     */
 /*                                                              */
 /*          omega   = proportion of NC items in group           */
 /*                                                              */
 /*          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.      */
 /*          (1990), Statistical Effects of Imperfect Inspection */
 /*          Sampling:  IV. Modified Dorfman Screening Pro-      */
 /*          cedures, Journal of Quality Technology 22, 128-137. */
 /*          See table on page 133.                              */
 /*                                                              */
 /*          Johnson, N. L., Kotz, S., and Wu, X. (1991).        */
 /*          Inspection Errors for Attributes in Quality         */
 /*          Control.  London:  Chapman & Hall.  See Chapter 9.  */
 /*                                                              */
 /****************************************************************/

data table;

   keep pcc pcnc omega k n0 n1 n2 n3 r10 r20 p pp outersum pct;

   format p 3.2 pp 3.2 pct 5.2 omega 4.2 pcc 5.4 pcnc 5.4 ;

   label outersum = 'Expected/Number/of Tests'
         pcc      = 'PC(C)'
         pcnc     = 'PC(NC)'
         k        = 'k'
         omega    = 'omega'
         p        = 'p'
         pp       = 'p'''
         pct      = 'Percent/Reduction'
         n0       = 'n0'
         n1       = 'n1'
         n2       = 'n2'
         n3       = 'n3'
         r10      = 'r10'
         r20      = 'r20'
         ;

   /*---set basic parameters---*/
   do omega=0.05 to 0.10 by 0.05;
   do k=0 to 2 by 1;
   do r10 = 1 to 2 by 1;
   do r20 = 1 to 2 by 1;

   r11 = 1; r21 = 1;
   r12 = 1; r22 = 1;
   r13 = 1; r23 = 1;

   /*---set parameters for group sizes---*/
   n0 = 12; h0 = 1;
   if k=0 then do;
      n1 = 1; h1 = 12;
      n2 = .; h2 = .;
      n3 = .; h3 = .;
      end;
   else
   if k=1 then do;
      n1 = 4; h1 = 3;
      n2 = 1; h2 = 4;
      n3 = .; h3 = .;
      end;
   else
   if k=2 then do;
      n1 = 4; h1 = 3;
      n2 = 2; h2 = 2;
      n3 = 1; h3 = 2;
   end;

   /*---loop over p error probabilities---*/
   do p=0.90 to 0.95 by 0.05;

      /*---compute p0* ---*/
      if (p<1.0) then do;

         p0s = probbnml(1-p,r10+r20-1,r20-1);
         p1s = probbnml(1-p,r11+r21-1,r21-1);
         p2s = probbnml(1-p,r12+r22-1,r22-1);
         p3s = probbnml(1-p,r13+r23-1,r23-1);
         end;

      else do;

         p0s = 1;
         p1s = 1;
         p2s = 1;
         p3s = 1;
         end;

      /*---loop over p' error probabilities---*/
      do pp=0.05 to 0.10 by 0.05;

         if pp>0.0 then do;
            p0ps = probbnml(1-pp,r10+r20-1,r20-1);
            p1ps = probbnml(1-pp,r11+r21-1,r21-1);
            p2ps = probbnml(1-pp,r12+r22-1,r22-1);
            p3ps = probbnml(1-pp,r13+r23-1,r23-1);
            end;
         else do;
            p0ps = 0;
            p1ps = 0;
            p2ps = 0;
            p3ps = 0;
            end;

         /*---sum over i ---*/
         outersum = 0.0;

         do i=0 to k+1;

            /*---evaluate product---*/
            if      i=0 then iproduct = 1;
            else if i=1 then iproduct = h1;
            else if i=2 then iproduct = h1*h2;
            else if i=3 then iproduct = h1*h2*h3;

            /*---evaluate first inner term---*/
            iterm1 = 0.0;
            do u=0 to i;

               /* evaluate diff = p0(nu) - p0(nu-1) */
               if      u=0 then diff=(1-omega)**n0;
               else if u=1 then diff=(1-omega)**n1 - (1-omega)**n0;
               else if u=2 then diff=(1-omega)**n2 - (1-omega)**n1;
               else if u=3 then diff=(1-omega)**n3 - (1-omega)**n2;
               else put 'error in diff for u > 3';

               /* evaluate prod1 = product of pv* */
               if      u=0 then prod1=1;
               else if u=1 then prod1=p0s;
               else if u=2 then prod1=p0s*p1s;
               else if u=3 then prod1=p0s*p1s*p2s;
               else put 'error in prod1 for u > 3';

               /*---evaluate prod2 = product of pw'* ---*/
               if i=0 then do;
                  if      u=0 then prod2 = 1.0 ;
                  else put 'error in prod2 for i=0';
                  end;
               else if i=1 then do;
                  if      u=0 then prod2 = p0ps;
                  else if u=1 then prod2 = 1;
                  else put 'error in prod2 for i=1';
                  end;
               else if i=2 then do;
                  if      u=0 then prod2 = p0ps*p1ps;
                  else if u=1 then prod2 = p1ps;
                  else if u=2 then prod2 = 1;
                  else put 'error in prod2 for i=2';
                  end;
               else if i=3 then do;
                  if      u=0 then prod2 = p0ps*p1ps*p2ps;
                  else if u=1 then prod2 = p1ps*p2ps;
                  else if u=2 then prod2 = p2ps;
                  else if u=3 then prod2 = 1;
                  else put 'error in prod2 for i=3';
                  end;
               else put 'error for i > 4';

               iterm1 = iterm1 + diff*prod1*prod2;
               end;   /* finish summing over u */

            /*---find eip = ei' ---*/
            if i=0 then do;
               if pp>0.0 then
                  eip = r20*probbnml(pp  ,r10+r20,r10-1)/(1-pp) +
                        r10*probbnml(1-pp,r10+r20,r20-1)/pp;
               else
                  eip = r20 ;
               end;
            else
            if i=1 then do;
               if pp>0.0 then
                  eip = r21*probbnml(pp  ,r11+r21,r11-1)/(1-pp) +
                        r11*probbnml(1-pp,r11+r21,r21-1)/pp;
               else
                  eip = r21 ;
               end;
            else if i=2 then do;
               if (pp>0.0) then
                  eip = r22*probbnml(pp  ,r12+r22,r12-1)/(1-pp) +
                        r12*probbnml(1-pp,r12+r22,r22-1)/pp;
               else
                  eip = r22 ;
               end;
            else
            if i=3 then do;
               if pp>0.0 then
                  eip = r23*probbnml(pp  ,r13+r23,r13-1)/(1-pp) +
                        r13*probbnml(1-pp,r13+r23,r23-1)/pp;
               else
                  eip = r23 ;
               end;
            else put 'error in eip for i > 3';

            iterm1 = iterm1*eip;

            /*---evaluate comp = 1-p0(ni) ---*/
            if      i=0 then comp = 1 - (1 - omega)**n0;
            else if i=1 then comp = 1 - (1 - omega)**n1;
            else if i=2 then comp = 1 - (1 - omega)**n2;
            else if i=3 then comp = 1 - (1 - omega)**n3;
            else put 'error in comp for i > 3';

            /*---evaluate prod3 = product of pv'* ---*/
            if      i=0 then prod3 = 1;
            else if i=1 then prod3 = p0s;
            else if i=2 then prod3 = p0s*p1s;
            else if i=3 then prod3 = p0s*p1s*p2s;
            else put 'error in prod3 for i > 3';

            /*---find ei = ei---*/
            if i=0 then do;
               if (p<1.0) then
                  ei = r20*probbnml(p  ,r10+r20,r10-1)/(1-p) +
                       r10*probbnml(1-p,r10+r20,r20-1)/p;
               else
                  ei = r10;
               end;
            else
            if i=1 then do;
               if (p<1.0) then
                  ei = r21*probbnml(p  ,r11+r21,r11-1)/(1-p) +
                       r11*probbnml(1-p,r11+r21,r21-1)/p;
               else
                  ei = r11;
               end;
            else
            if i=2 then do;
               if (p<1.0) then
                  ei = r22*probbnml(p  ,r12+r22,r12-1)/(1-p) +
                       r12*probbnml(1-p,r12+r22,r22-1)/p;
               else
                  ei = r12;
               end;
            else
            if i=3 then do;
               if (p<1.0) then
                  ei = r23*probbnml(p  ,r13+r23,r13-1)/(1-p) +
                       r13*probbnml(1-p,r13+r23,r23-1)/p;
               else
                  ei = r13;
               end;
            else put 'error in prod3 for i > 3';

            iterm2 = comp*prod3*ei;

            /*---update outer sum---*/
            outersum = outersum + iproduct*(iterm1 + iterm2);

            end;  /* over i */

         /*---percent savings---*/
         pct = 100.0*(1 - outersum/n0);

         /* pcnc = pc( nc ) */
         if k = 0 then pcnc = p0s * p1s;
         else
         if k = 1 then pcnc = p0s * p1s * p2s;
         else
         if k = 2 then pcnc = p0s * p1s * p2s * p3s;
         else put 'error for pcnc with k > 2';

         /*---pcc = pc( c )---*/
         sumpcc = 0.0;
         do s = 0 to k + 1;

            /*---compute diff = p0*(ns) - p0*(ns-1)---*/
            if s = 0 then diff = (1-omega)**(n0-1);
            else
            if s = 1 then do;
               if s=k+1 then
                  diff = 1                -(1-omega)**(n0-1);
               else
                  diff = (1-omega)**(n1-1)-(1-omega)**(n0-1);
               end;
            else
            if s = 2 then do;
               if s=k+1 then
                  diff = 1                -(1-omega)**(n1-1);
               else
                  diff = (1-omega)**(n2-1)-(1-omega)**(n1-1);
               end;
            else
            if s = 3 then do;
               if s=k+1 then
                  diff = 1                -(1-omega)**(n2-1);
               else
                  diff = (1-omega)**(n3-1)-(1-omega)**(n2-1);
               end;
            else put 'error for pcc diff with s > 3';

            /*---prod1 = product of pj* ---*/
            if      s = 0 then prod1 = 1;
            else if s = 1 then prod1 = p0s ;
            else if s = 2 then prod1 = p0s * p1s ;
            else if s = 3 then prod1 = p0s * p1s * p2s ;
            else put 'error for pcc prod1 with s > 3 ';

            /*---prod2 = product of pj* ---*/
            if      s = 0 then do;
               if      k = 0 then prod2 = p0ps * p1ps;
               else if k = 1 then prod2 = p0ps * p1ps * p2ps;
               else if k = 2 then prod2 = p0ps * p1ps * p2ps *p3ps;
               else put 'error with s = 0 and k > 2';
               end;
            else if s = 1 then do;
               if      k = 0 then prod2 = p1ps;
               else if k = 1 then prod2 = p1ps * p2ps;
               else if k = 2 then prod2 = p1ps * p2ps * p3ps;
               else put 'error with s = 1 and k > 2';
               end;
            else if s = 2 then do;
               if      k = 0 then prod2 = 1;
               else if k = 1 then prod2 = p2ps;
               else if k = 2 then prod2 = p2ps * p3ps;
               else put 'error with s = 2 and k > 2';
               end;
            else if s = 3 then do;
               if      k = 0 then prod2 = 1;
               else if k = 1 then prod2 = 1;
               else if k = 2 then prod2 = p3ps;
               else put 'error with s = 3 and k > 2';
               end;
            else put 'error for pcc prod2 with s > 3 ';

            sumpcc = sumpcc + diff * prod1 * prod2;
            end; /* over s */

         pcc = 1 - sumpcc;

         output;

         end;   /* over pp */
      end;   /* over p */

   end;  /* over r10 */
   end;  /* over r20 */
   end;  /* over k */
   end;  /* over omega */
   run;

proc sort; by n0 omega k;

proc print label noobs split='/' ;
   var n1 n2 n3 r10 r20 p pp pct pcc pcnc;
   by n0 omega k;
run;