#
# Written by:
# -- 
# John L. Weatherwax                2009-04-21
# 
# email: wax@alum.mit.edu
# 
# Please send comments and especially bug reports to the
# above email address.
#
#-----

source('utils.R')

# Ex 16.31:
#
c = 2
n = 50
N = 500
xs = 0:c
ps = seq( 0.01, 0.1, by=0.01 )
P_A_hyper = rep( NA, length(ps) ) # use the hypergeometric distribution
P_A_binom = rep( NA, length(ps) ) # use the binomial distribution
for( pi in 1:length(ps) ){
    p = ps[pi]
    P_A_hyper[pi] = sum( dhyper( xs, round(N*p), round(N*(1-p)), n ) )
    P_A_binom[pi] = sum( dbinom( xs, n, p ) )
}

#postscript("../../WriteUp/Graphics/Chapter16/ex_31.eps", onefile=FALSE, horizontal=FALSE)
plot( ps, P_A_hyper, type='l', col='green', pch=19, xlab='p', ylab='P(A)' )
lines( ps, P_A_binom, type='l', col='red' )
grid()
legend( 0.07, 1.0, c('hypergeometric','binomial'), lty=c(1, 1), col=c('green', 'red') )
#dev.off()



# Ex 16.32:
# 
n = 50
N = 5000
c = 1
xs = 0:c
ps = seq( 0.01, 0.1, by=0.005 )
P_A = rep( NA, length(ps) )
for( pi in 1:length(ps) ){
    p = ps[pi]
    P_A[pi] = sum( dhyper( xs, round(N*p), round(N*(1-p)), n ) )
}
#postscript("../../WriteUp/Graphics/Chapter16/ex_32.eps", onefile=FALSE, horizontal=FALSE)
plot( ps, P_A, type='l', col='green', pch=19, xlab='p', ylab='P(A)' )
grid()
#dev.off()


# Ex 16.33:
n = 100
c = 2
xs = 0:c

P2_A = rep( NA, length(ps) ) # a second possible plan
for( pi in 1:length(ps) ){
    p = ps[pi]
    P2_A[pi] = sum( dhyper( xs, round(N*p), round(N*(1-p)), n ) )
}

#postscript("../../WriteUp/Graphics/Chapter16/ex_32.eps", onefile=FALSE, horizontal=FALSE)
plot( ps, P_A, type='l', col='green', pch=19, xlab='p', ylab='P(A)' )
lines( ps, P2_A, type='l', col='red', pch=19, xlab='p', ylab='P(A)' )
grid()
legend( 0.07, 0.92, c('n=50; c=1','n=100; c=2'), lty=c(1, 1), col=c('green', 'red') )
#dev.off()


# 16.34:
#
AQL = 0.02 # p_1
LTPD = 0.07 # p_2
LTPD / AQL
c = 5
n = 2.613 / AQL
n = ceiling(n)

xs = 0:c
alpha = 1 - sum( dbinom( xs, n, 0.02 ) )
beta = sum( dbinom( xs, n, 0.07 ) )
print( c(alpha, beta) )



# 16.35:
#
n1 = 50
n2 = 50
c1 = 1
r1 = 4
c2 = 5
PAs = c()
for( p in c( 0.01, 0.05, 0.1 ) ){
    P_A = sum( dbinom( 0:01, n1, p ) ) +
          dbinom( 2, n1, p ) * sum( dbinom( 0:3, n2, p ) ) +
          dbinom( 3, n1, p ) * sum( dbinom( 0:2, n2, p ) )
}
PAs = c( PAs, P_A )



# 16.36:
#
n1 - 50
n2 = 100
c1 = 1
r1 = 4
c2 = 3
PAs = c()
for( p in c( 0.02, 0.05, 0.1 ) ){
    P_A = sum( dbinom( 0:c1, n1, p ) ) +
          dbinom( 2, n1, p ) * sum( dbinom( 0:1, n2, p ) ) +
          dbinom( 3, n1, p ) * dbinom( 0, n2, p )
    PAs = c( PAs, P_A )
}
print( PAs ) 



# 16.37:
#
n1 = 50
c = 2
ps = seq( 0.01, 0.1, by=0.01 )
AOQ = ps * pbinom( 0, n1, ps )

N = 2000
ATI = n1 * pbinom( c, n1, ps ) + N * ( 1 - pbinom( c, n1, ps ) )


# 16.38:
#
n1 = 50
c = 1
ps = seq( 0.01, 0.1, by=0.01 )
AOQ = ps * pbinom( c, n1, ps )

N = 2000
ATI = n1 * pbinom( c, n1, ps ) + N * ( 1 - pbinom( c, n1, ps ) )
print( max( AOQ ) )