PC(C) for Curtailed Group Testing-Finite
/****************************************************************/
/* S A S S A M P L E L I B R A R Y */
/* */
/* NAME: IECURT3 */
/* TITLE: PC(C) for Curtailed Group Testing-Finite */
/* PRODUCT: QC */
/* SYSTEM: ALL */
/* KEYS: Inspection Sampling, */
/* PROCS: TABULATE */
/* DATA: */
/* */
/* MISC: */
/* */
/* NOTES: This program tabulates the probability of correct */
/* classification of a conforming item 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.3. */
/* */
/* 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 = 'PC(C)';
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 - 1, d );
lognum = lgamma( d + 1 )
+ lgamma( nlot - d )
+ lgamma( nsample )
+ lgamma( ( nlot - nsample ) + 1 );
lgnlot = lgamma( nlot );
/*---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 = 1 to nsample;
do z2 = 0 to nsample - 1;
/*---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 - 1;
link cond;
fact2 = cprob;
/*---conditional probability---*/
fact3 = ( 1 - pprime ) *
min( ( nsample-z1 )/( nsample-z2 ), 1 );
fact3 = fact3 +
pprime * max((1+z2-z1) / (z2+1),
0 );
innersum = innersum + fact1 * fact2 * fact3;
end;
end;
/*---hypergeometric probability---*/
logden = lgamma( y + 1 )
+ lgamma( ( d - y ) + 1 )
+ lgamma( nsample - y )
+ lgamma( ( nlot-d ) - ( nsample-y ) + 1 )
+ lgnlot;
hypprob = exp( lognum - logden );
outersum = outersum + innersum * hypprob;
end; /* finish loop over y values */
result = 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;
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 ( q>0 ) & ( 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;
/*---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;
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 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;