# 
# 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.
# 
#-----

N = 1
N = 10
pr_a_S1 = 0.55
pr_b_S1 = 1-pr_a_S1
pr_a_S2 = 0.58
pr_b_S2 = 1-pr_a_S2
  
gamma_I = seq( from=0, to=N+1 ) 

P_F = 1 - pbinom( gamma_I-1, N, pr_a_S1 )
P_D = 1 - pbinom( gamma_I-1, N, pr_a_S2 )

#postscript("../../WriteUp/Graphics/Chapter2/prob_2.2.7_N_1.eps", onefile=FALSE, horizontal=FALSE)
#postscript("../../WriteUp/Graphics/Chapter2/prob_2.2.7_N_10.eps", onefile=FALSE, horizontal=FALSE)

plot( P_F, P_D, col="black", type="p", xlim=c(0,1), ylim=c(0,1) )
lines( gamma, gamma, col="gray", type="l" )
abline( v=0.25, col="gray" )
#dev.off()

# Solve for p_gamma_2
#

# Write P_F (and all other variables) in increasing order:
#
pf_order = order(P_F)
P_F = P_F[pf_order]
P_D = P_D[pf_order]
gamma_I = gamma_I[pf_order]

alpha = 0.25

res = hist( alpha, breaks=P_F, plot=FALSE )

bin = which( res$counts==1 )
print( bin )
print( gamma_I[bin] ) 

p_gamma_2 = ( alpha - P_F[bin+1] )/(P_F[bin]-P_F[bin+1])
print( p_gamma_2 )