E(N) for Curtailed Group Testing-Infinite
/****************************************************************/
/* S A S S A M P L E L I B R A R Y */
/* */
/* NAME: IECURT1B */
/* TITLE: E(N) for Curtailed Group Testing-Infinite */
/* 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 an infinite lot */
/* under an imperfect inspection model. */
/* */
/* Notation: */
/* */
/* omega = proportion of defective 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 10.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 omega nsample p1 p1prime p pprime result;
label omega = 'omega'
nsample = 'n (sample)'
p1 = 'p1'
p1prime = 'p1'''
p = 'p'
pprime = 'p'''
result = 'E(Number Tests)';
format result 6.4 ;
/*---set main parameters---*/
omega = 0.05;
nsample = 5;
p1 = 0.95;
p1prime = 0.05;
/*---local variables---*/
ylower = 0;
yupper = nsample;
/*---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;
nval = nsample;
link cond;
fact1 = cprob;
*--second conditional probability factor--;
q = p;
qprime = pprime;
zval = z2;
yval = y;
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;
/*---binomial probability---*/
p_ = omega;
n_ = nsample;
k_ = y;
link binomial;
outersum = outersum + innersum * binprob;
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 */
/* */
/* nval = sample 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;
if ( q = 0 ) & ( qprime = 0 ) then do;
if zval = 0 then cprob = 1;
end;
else
if ( q = 0 ) & ( abs( qprime - 1 ) < fuzz ) then do;
if zval = nval - yval then cprob = 1;
end;
else
if ( abs( q - 1 ) < fuzz ) & ( qprime = 0 ) then do;
if zval = yval then cprob = 1;
end;
else
if ( abs( q - 1 ) < fuzz ) & ( qprime > 0 ) & ( qprime < 1 )
then do;
n_ = nval - yval;
p_ = qprime;
k_ = zval - yval;
link binomial;
cprob = binprob;
end;
else
if ( 0 < q ) & ( q < 1 ) & ( abs( qprime - 1 ) < fuzz ) then do;
n_ = yval;
p_ = q;
k_ = zval - ( nval - yval );
link binomial;
cprob = binprob;
end;
else
if ( abs(q-1) < fuzz ) & ( abs(qprime-1) < fuzz ) then do;
if zval =nval then cprob = 1;
end;
else
if ( 0 < q ) & ( q < 1) & ( qprime > 0 ) & ( qprime < 1 )
then do;
if ( 0 <= zval ) & ( zval <= nval ) then do;
/*---convolution of binomial distributions---*/
do xval = 0 to zval 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;
else
if ( q = 0 ) & ( 0 < qprime ) & ( qprime < 1 ) then do;
n_ = nval - yval;
p_ = qprime;
k_ = zval;
link binomial;
cprob=binprob;
end;
else
if ( qprime = 0 ) & ( 0 < q ) & ( q < 1 ) then do;
n_ = yval;
p_ = q;
k_ = zval;
link binomial;
cprob = binprob ;
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;
run;
proc sort data=table;
by omega nsample p1 p1prime;
proc tabulate data=table noseps;
by omega nsample p1 p1prime;
class p pprime;
var result;
table p, pprime*result=' '*sum=' '*f=8.4 / rts=7;
run;