2PC(C)n|y for Two-Stage Dorfman-Sterrett
/****************************************************************/
/* S A S S A M P L E L I B R A R Y */
/* */
/* NAME: IEDS2C */
/* TITLE: 2PC(C)n|y for Two-Stage Dorfman-Sterrett */
/* PRODUCT: QC */
/* SYSTEM: ALL */
/* KEYS: Inspection Sampling, */
/* PROCS: TABULATE */
/* DATA: */
/* */
/* MISC: */
/* */
/* NOTES: This program tabulates the conditional probability */
/* of correct classification 2PC(C)n|y for the two- */
/* stage Dorfman-Sterrett procedure under an imperfect */
/* inspection model. */
/* */
/* Notation: */
/* */
/* n = number of items in group */
/* y = number of truly NC items in group */
/* */
/* p0 = Pr[ for group test, declare NC | at least */
/* one NC item ] */
/* */
/* p0prime = 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 ] */
/* */
/* Note that p0prime does not enter into the */
/* computations. */
/* */
/* */
/* REF: Johnson, N. L., Kotz, S., and Rodriguez, R. N. */
/* (1990), Statistical Effects of Imperfect Inspection */
/* Sampling: IV. Modified Dorfman Screening Pro- */
/* cedures, Journal of Quality Technology 22, 128-137. */
/* See Table 11.3. */
/* */
/* Johnson, N. L., Kotz, S., and Wu, X. (1991). */
/* Inspection Errors for Attributes in Quality */
/* Control. London: Chapman & Hall. See Chapter 7. */
/* */
/****************************************************************/
data table;
keep k _y_ nsample p0 p0prime p pprime result;
label nsample = 'n (sample)'
k = 'k (stage)'
p0 = 'p0'
p0prime = 'p0'''
p = 'p'
pprime = 'p'''
result = '2PC(C)n|y'
_y_ = 'y' ;
/*---set main parameters---*/
k = 2;
nsample = 10;
p0 = 0.75;
p0prime = 0.05;
/*---loop over y values---*/
do _y_ = 0 to nsample - 1;
/*---loop over p values---*/
do p = 0.75, 0.95;
/*---loop over pprime values---*/
do pprime = 0.05, 0.25;
/*---special handling for _y_=0---*/
if _y_ = 0 then do;
outsum = 0.0;
upper_m_ = nsample-2;
do _m_=1 to upper_m_;
/*---find probpcc1 = 1pc(c)n-m|0---*/
n = nsample-_m_;
y = 0;
link prbccfun;
outsum = outsum
+ (1 + (nsample-_m_)*(1-probpcc1))*pprime*
(1-pprime)**(_m_-1);
end;
probpcc2 = outsum + 2*pprime*(1-p)**(nsample-2);
probpcc2 = 1 - p0prime*probpcc2/nsample;
end;
/*---general _y_---*/
else do;
/*---set upper limit for m---*/
upper_m_ = nsample - 2;
/*---find outer sum---*/
outsum = 0.0;
do _m_=1 to upper_m_ by 1;
/*---find inner sum---*/
insum = 0.0;
lower_t_ = max( 0, _y_+_m_-nsample);
upper_t_ = min( _m_, _y_);
do _t_=lower_t_ to upper_t_ by 1;
/*---find pc2 = pc(m,t|n,y)---*/
m = _m_;
t = _t_;
n = nsample;
y = _y_;
link pcfun;
pc2 = pc;
/*---find pnc = pnc(m,t|n,y)---*/
link pncfun;
/*---find pmtny2 = p(m,t|n,y)---*/
pmtny2 = pnc + pc;
/*---find probpcc1 = 1pc(c)n-m|y-t---*/
n = nsample-_m_;
y = _y_-_t_;
link prbccfun;
/*---increment inner sum---*/
insum = insum
+ pc2
+ pmtny2 *
( nsample + _t_ - ( _m_+_y_ ) ) *
( 1 - probpcc1 );
end; /* finish loop for inner sum */
/*---increment outer sum---*/
outsum = outsum + insum;
end; /* finish loop for outer sum */
/*---find _term2_ = pc(n-1,y|n,y)---*/
m = nsample-1;
t = _y_;
n = nsample;
y = _y_;
link pcfun;
_term2_ = pc;
/*---find _term3_ = p'p(n-1,y|n,y)---*/
link pncfun;
_term3_ = pprime*(pc + pnc);
/*---find _term4_ = pc(n,y|n,y)---*/
m = nsample;
t = _y_;
n = nsample;
y = _y_;
link pcfun;
_term4_ = pc;
/*---find _term5_ = pc(n-1,y-1|n,y)---*/
m = nsample-1;
t = _y_-1;
n = nsample;
y = _y_;
link pcfun;
_term5_ = pc;
probpcc2 = 1 -
p0*(outsum+_term2_+_term3_+_term4_+_term5_)/
(nsample-_y_);
end; /* finish for general y */
result = probpcc2;
output;
end; /* finish loop over pprime values */
end; /* finish loop over p values */
end; /* finish loop over _y_ values */
return; /* finish main program */
/*------------------------------------------------------------*/
/* */
/* Calculate 1PC(C)n|y for y = 0, ..., n - 1 */
/* */
/*------------------------------------------------------------*/
prbccfun:
/*---initialize---*/
probpcc1 = 0.0 ;
/*---special handling if y=0---*/
if y = 0 then do;
probpcc1 = 1 + n*pprime*p0prime - p0prime;
probpcc1 = probpcc1
- (1-p0prime)*(1-2*pprime)*(1-pprime)**(n-2);
probpcc1 = 1 - p0prime*probpcc1/n;
end;
else if y = n then probpcc1 = 0.0 ;
/*---general case: y>0---*/
else do;
/*---set upper limit for m---*/
upperm=n-2;
/*---find outer sum---*/
outersum=0.0;
do m=1 to upperm by 1;
/*---find inner sum---*/
innersum=0.0;
lowert=max(0,y+m-n);
uppert=min(m,y-1);
do t=lowert to uppert by 1;
/*---find pc = pc(m,t|n,y)---*/
link pcfun;
/*---find pnc = pnc(m,t|n,y)---*/
link pncfun;
/*---find pmtny = p(m,t|n,y)---*/
pmtny = pnc + pc;
/*---increment inner sum---*/
innersum = innersum
+ (pc + pmtny*((n+t)-(m+y))*p0*pprime);
end; /* finish loop for inner sum */
/*---find pc(m,y|n,y)---*/
t=y;
link pcfun;
/*---find pnc(m,y|n,y)---*/
t=y;
link pncfun;
/*---increment outer sum---*/
outersum = outersum
+ innersum
+ pc
+ (pc + pnc)*(n-m)*p0prime*pprime;
end; /* finish loop for outer sum */
/*---find term2 = pc(n-1,y|n,y)---*/
m=n-1;
t=y;
link pcfun;
term2 = pc;
/*---find term3 = p'(pnc(n-1,y|n,y) + pc(n-1,y|n,y))---*/
m=n-1;
t=y;
link pncfun;
term3 = pprime*(pnc+pc);
/*---find term4 = pc(n,y|n,y)---*/
m=n;
t=y;
link pcfun;
term4 = pc;
/*---find term5 = pc(n-1,y-1|n,y)---*/
m=n-1;
t=y-1;
link pcfun;
term5 = pc;
probpcc1=1 -
p0*(outersum + term2 + term3 + term4 + term5)/(n-y);
end;
return; /* finish computing 1PC(C)n|y */
/*------------------------------------------------------------*/
/* */
/* Calculate pnc = PNC( m, t | n, y ) */
/* */
/*------------------------------------------------------------*/
pncfun:
if t<=0 then pnc = 0.0;
else if t>m then pnc = 0.0;
else do;
if p = 1.0 then do;
if t = 1 then do;
lognum = lgamma( 1 + n - m )
+ ( m - 1 ) * log( 1 - pprime )
+ lgamma( y + 1 )
+ lgamma( 1 + n - y );
logden = lgamma( y )
+ lgamma( 1 + ( n - m ) - ( y - 1 ) )
+ lgamma( 1 + n );
pnc = exp( lognum - logden );
end;
else pnc = 0.0;
end;
else do;
lognum = lgamma( m )
+ lgamma( 1 + n - m )
+ ( t - 1 ) * log( 1 - p )
+ log( p )
+ ( m - t ) * log( 1 - pprime )
+ lgamma( y + 1 )
+ lgamma( 1 + n - y );
logden = lgamma( t )
+ lgamma( 1 + m - t )
+ lgamma( 1 + y - t )
+ lgamma( 1 + ( n - m ) - ( y - t ) )
+ lgamma( 1 + n );
pnc = exp( lognum - logden );
end;
end;
return; /* finish calculating pnc */
/*------------------------------------------------------------*/
/* */
/* Calculate pc = PC( m, t | n, y ) */
/* */
/*------------------------------------------------------------*/
pcfun:
if t >= m then pc = 0.0;
else if t < 0 then pc = 0.0;
else do;
if pprime = 0 then pc = 0.0;
else if p = 1 then do;
if t = 0 then do;
lognum = lgamma( 1 + n - m )
+ ( m - 1 ) * log( 1 - pprime )
+ log( pprime )
+ lgamma( y + 1 )
+ lgamma( 1 + n - y );
logden = lgamma( 1 + y )
+ lgamma( 1 + ( n - m ) - y )
+ lgamma( 1 + n );
pc = exp( lognum - logden );
end;
else pc = 0.0;
end;
else do;
lognum = lgamma( m )
+ lgamma( 1 + n - m )
+ t * log( 1 - p )
+ (m - t - 1 ) * log( 1 - pprime )
+ log( pprime )
+ lgamma( y + 1 )
+ lgamma( 1 + n - y );
logden = lgamma( 1 + t )
+ lgamma( m - t )
+ lgamma( 1 + y - t )
+ lgamma( 1 + ( n - m ) - ( y - t ) )
+ lgamma( 1 + n );
pc = exp( lognum - logden );
end;
end;
return; /* finish calculating pc */
run;
proc sort data=table;
by k nsample p0 p0prime p _y_ pprime;
proc tabulate data=table noseps;
by k nsample p0 p0prime p;
class _y_ pprime;
var result;
table _y_, pprime*result=' '*sum=' '*f=8.4 / rts=5;
run;
