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;