# Bayesian sequencial learning 
#
# theta=1 if FLU and theta=0 otherwise.
#
# Not-so-good clinical test X
# Pr(x=1|theta=1) = 0.9
# Pr(x=1|theta=0) = 0.7
#
# Much-better clinical test Y 
# Pr(y=1|theta=1) = 0.9
# Pr(y=1|theta=0) = 0.1
#
# Prevalence of disease
# Pr(theta=1) = alpha
#
# Therefore, 
# 
# Pr(theta=1|x=1) =  1/(1+(1-alpha)/alpha*1/LR)
#
# where LR = Pr(x=1|theta=1)/Pr(x=1|theta=0) is the likelihood ratio.

px1t1=0.9
px1t0=0.7
py1t1=0.9
py1t0=0.1
alpha = 0.1

LR.x1 = px1t1/px1t0
LR.y1 = py1t1/py1t0
LR.x1y1 =  LR.x1*LR.y1

prob.x1 = 1/(1+(1-alpha)/alpha*(1/LR.x1))
prob.y1 = 1/(1+(1-alpha)/alpha*(1/LR.y1))
prob.x1y1 = 1/(1+(1-alpha)/alpha*(1/LR.x1y1))
c(prob.x1,prob.y1,prob.x1y1)

par(mfrow=c(1,3))
for (px1t0 in c(0.7,0.4,0.2)){
LR.x1 = px1t1/px1t0
LR.y1 = py1t1/py1t0
LR.x1y1 =  LR.x1*LR.y1

alpha = seq(0,1,length=100)
plot(alpha,1/(1+(1-alpha)/alpha*(1/LR.x1)),col=2,type="l",xlab="Prior of theta=1",ylab="Posterior of theta=1")
lines(alpha,1/(1+(1-alpha)/alpha*(1/LR.y1)),col=3)
lines(alpha,1/(1+(1-alpha)/alpha*(1/LR.x1y1)),col=4)
abline(0,1)
title(paste("Pr(x=1|theta=0)=",px1t0,sep=""))
}
legend("bottomright",legend=c("Pr(theta=1|x=1)","Pr(theta=1|y=1)","Pr(theta=1|x=1,y=1)"),col=c(2,3,4),lty=1)