Acceptance Probabilities for Partial Link Sampling

``` /****************************************************************/
/*          S A S   S A M P L E   L I B R A R Y                 */
/*                                                              */
/*   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---*/

/*---computer term2---*/

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;

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;
product = zprob;

/* obtain joint probability that Z2=z2, Z2P=z2p */
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;
hprob1 = hypprob;

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

mhyp = hprob1 * hprob2;

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

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

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;

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

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

p_ = pprime;
k_ = zval - xval;
n_ = nval - yval;
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;

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;

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;

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;

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

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;

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;

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

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;

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

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

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

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

/*---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;
factor1 = binprob;

p_ = pprime;
k_ = z - x;
n_ = nsample - y;
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;

```