E(N) for Curtailed Group Testing-Finite
/****************************************************************/
/* S A S S A M P L E L I B R A R Y */
/* */
/* NAME: IECURT1 */
/* TITLE: E(N) for Curtailed Group Testing-Finite */
/* PRODUCT: QC */
/* SYSTEM: ALL */
/* KEYS: Inspection Sampling, */
/* PROCS: TABULATE */
/* DATA: */
/* */
/* MISC: */
/* */
/* NOTES: This program tabulates the expected number of tests */
/* for curtailed group testing from a finite lot under */
/* an imperfect inspection model. */
/* */
/* Notation: */
/* */
/* nlot = size of lot */
/* d = number of nonconforming items in lot */
/* */
/* n = sample size */
/* */
/* p1 = Pr[ for group test, declare NC | at least */
/* one NC item ] */
/* */
/* p1prime = Pr[ for group test, declare NC | no NC */
/* items ] */
/* */
/* 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. */
/* (1986), Statistical Effects of Imperfect Inspection */
/* Sampling: III. Screening (Group Testing), Journal */
/* of Quality Technology 20, 98-124. See Table 9.1. */
/* */
/* Johnson, N. L., Kotz, S., and Wu, X. (1991). */
/* Inspection Errors for Attributes in Quality */
/* Control. London: Chapman & Hall. See Chapter 8. */
/* */
/****************************************************************/
data table;
keep nlot d nsample p1 p1prime p pprime result;
label nlot = 'N (lot)'
d = 'D'
nsample = 'n (sample)'
p1 = 'p1'
p1prime = 'p1'''
p = 'p'
pprime = 'p'''
result = 'E(Number Tests)';
format result 6.4 ;
/*---set main parameters---*/
nlot = 100;
nsample = 5;
d = 5;
p1 = 0.95;
p1prime = 0.05;
/*---local variables---*/
ylower = max( 0, nsample + d - nlot );
yupper = min( nsample, d );
/*---loop over p values---*/
do p = 0.75, 0.80, 0.85, 0.90, 0.95;
/*---loop over pprime values---*/
do pprime = 0.0, 0.025, 0.05, 0.075, 0.10 ;
/*---outer sum---*/
outersum = 0.0 ;
do y = ylower to yupper by 1;
/*---inner double sum---*/
innersum = 0.0;
do z1 = 0 to nsample;
do z2 = 0 to nsample;
/*---first conditional probability factor---*/
q = p1;
qprime = p1prime;
zval = z1;
yval = y;
dval = d;
nval = nsample;
link cond;
fact1 = cprob;
*--second conditional probability factor--;
q = p;
qprime = pprime;
zval = z2;
yval = y;
dval = d;
nval = nsample;
link cond;
fact2 = cprob;
/*---conditional expectation factor--*/
if z2 > z1 then
fact3 = (( nsample + 1 )*z1)/(z2+1);
else if z1 = z2 then
fact3 = z1*(nsample-z1)*((1/(z1+1)) +
(1/(1+nsample-z1)));
else
fact3=((nsample+1)*(nsample-z1))/
(nsample+1-z2);
innersum = innersum + fact1 * fact2 * fact3;
end;
end;
/*---hypergeometric probability---*/
bign_ = nlot;
litn_ = nsample;
d_ = d;
y_ = y;
link hypergmt;
outersum = outersum + innersum * hypprob;
end; /* finish loop over y values */
result = 1 + outersum ;
output;
end; /* finish loop over pprime values */
end; /* finish loop over p values */
return; /* finish main program */
/*------------------------------------------------------------*/
/* */
/* 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 */
/* q = Pr[ correctly classifying a defective item ] */
/* qprime = Pr[ incorrectly classifying a good item ] */
/* */
/* */
/*------------------------------------------------------------*/
cond:
/*---initialize result to zero---*/
cprob = 0.0;
fuzz = 0.001;
/*---set limits for subscript---*/
lolim = max( 0, nval + dval - nlot );
uplim = min( nval, dval );
if ( q = 0 ) & ( qprime = 0 ) then do;
if ( yval >= lolim ) & ( yval <= uplim) then do;
if zval = 0 then cprob = 1;
end;
end;
else
if ( q = 0 ) & ( abs( qprime - 1 ) < fuzz ) then do;
if ( yval >= lolim ) & ( yval <= uplim ) then do;
if zval = nval - yval then cprob = 1;
end;
end;
else
if ( abs( q - 1 ) < fuzz ) & ( qprime = 0 ) then do;
if ( lolim <= yval ) & ( yval <= uplim ) then do;
if zval = yval then cprob = 1;
end;
end;
else
if ( abs( q - 1 ) < fuzz ) & ( qprime > 0 ) & ( qprime < 1 )
then do;
if ( lolim <= yval ) & ( yval <= uplim ) then do;
if ( yval <= zval ) & ( zval <= nval ) then do;
n_ = nval - yval;
p_ = qprime;
k_ = zval - yval;
link binomial;
cprob = binprob;
end;
end;
end;
else
if ( q>0 ) & ( q<1 ) & ( abs( qprime - 1 ) < fuzz ) then do;
if ( lolim <= yval ) & ( yval <= uplim ) then do;
if ( nval - yval <= zval ) & ( zval <= nval ) then do;
n_ = yval;
p_ = q;
k_ = zval - ( nval - yval );
link binomial;
cprob = binprob;
end;
end;
end;
else
if ( abs(q-1) < fuzz ) & ( abs(qprime-1) < fuzz ) then do;
if ( lolim <= yval ) & ( yval <= uplim ) then do;
if zval =nval then cprob = 1;
end;
end;
else
if ( 0 < q ) & ( q < 1) & ( qprime > 0 ) & ( qprime < 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_ = q;
k_ = xval;
n_ = yval;
link binomial;
f1 = binprob;
p_ = qprime;
k_ = zval - xval;
n_ = nval - yval;
link binomial;
f2 = binprob;
cprob =cprob + f1 * f2;
end;
end;
end;
end;
else
if ( q = 0 ) & ( 0 < qprime ) & ( qprime < 1 ) then do;
if ( yval >= lolim ) & ( yval <= uplim ) then do;
if ( zval >= 0 ) & ( zval <= nval - yval ) then do;
n_ = nval - yval;
p_ = qprime;
k_ = zval;
link binomial;
cprob=binprob;
end;
end;
end;
else
if ( qprime = 0 ) & ( 0 < q ) & ( q < 1 ) then do;
if ( lolim <= yval ) & ( yval <= uplim ) then do;
if ( zval >= 0 ) & ( zval <= yval ) then do;
n_ = yval;
p_ = q;
k_ = zval;
link binomial;
cprob = binprob ;
end;
end;
end;
return; /* finish cond */
/*---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 p1 p1prime;
proc tabulate data=table noseps;
by nlot d nsample p1 p1prime;
class p pprime;
var result;
table p, pprime*result=' '*sum=' '*f=8.4 / rts=7;
run;