Distribution of Defective Items-Finite
/****************************************************************/
/* S A S S A M P L E L I B R A R Y */
/* */
/* NAME: IEDISTZ1 */
/* TITLE: Distribution of Defective Items-Finite */
/* PRODUCT: QC */
/* SYSTEM: ALL */
/* KEYS: Inspection Sampling, */
/* PROCS: TABULATE */
/* DATA: */
/* */
/* MISC: */
/* */
/* NOTES: This program tabulates the distribution of Z, the */
/* apparent number of nonconforming (defective) items */
/* in a sample drawn from a single, finite lot under */
/* an imperfect inspection model. */
/* */
/* Notation: */
/* */
/* zprob = Pr[ Z = z ] */
/* */
/* z = apparent number of nonconforming items */
/* nlot = size of lot */
/* d = number of nonconforming items in lot */
/* nsample = sample size */
/* 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. */
/* (1985), Statistical Effects of Imperfect Inspection */
/* Sampling: I. Some Basic Distributions, Journal of */
/* Quality Technology 17, 1-31. See Table 1. */
/* */
/* Johnson, N. L., Kotz, S., and Wu, X. (1991). */
/* Inspection Errors for Attributes in Quality */
/* Control. London: Chapman & Hall. See Chapter 2. */
/* */
/****************************************************************/
data table;
keep nlot d nsample p pprime z zprob;
label nlot = 'N (lot)'
nsample = 'n (sample)'
d = 'D'
p = 'p'
pprime = 'p'''
z = 'z'
zprob = 'Pr[ Z = z ]' ;
format zprob 6.4 ;
/*---set main parameters---*/
nlot = 100;
d = 5;
nsample = 5;
/*---used for roundoff---*/
fuzz = 0.0001 ;
/*---lower and upper limits for y---*/
miny = max( 0, nsample + d - nlot );
maxy = min( nsample, d );
/*---loop over range of p values---*/
do p = 0.75 to 1.00 by 0.05 ;
/*---loop over range of pprime values---*/
do pprime = 0.0 to 0.10 by 0.025 ;
/*---determine Pr[ Z = z ]---*/
do z = 0 to nsample by 1 ;
/*---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 */
/*---output Pr[ Z = z ] ---*/
output;
end; /* finish loop over z */
end; /* finish loop over pprime */
end; /* finish loop over p */
/*---finish main program---*/
return;
/*---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 p;
proc tabulate data=table noseps;
by nlot d nsample p;
class z pprime;
var zprob;
table z, pprime*zprob=' '*sum=' '*f=8.4 / rts=5;
run;