# Bayes factor example, discussed in class on Tue Feb 16
# H1: Y|theta ~ Binomial(20, theta) where theta = 0.5
# H2: Y|theta ~ Binomial(20, theta) where theta ~ Uniform(0,1)
rm(list=ls())
# First let’s do classical hypothesis test of
# Y ~ Binomial(n=20, theta)
# H0: theta = .5 versus H1: theta ^= .5
# for comparison
# For what y-values is p-value < .05? < .01?
y <- 0:20
p.value <- c(2*pbinom(0:9, 20, .5), 1, 2*(1-pbinom(10:19, 20, .5)))
plot(log10(p.value) ~ y, type="b")
abline(h=log10(.05), col="blue")
abline(h=log10(.01), col="red")
legend("bottom", inset=.10, lty=1, col=c("blue","red"),
legend=c("p-value<.05: y<=5 or y>=15″,
“p-value<.01: y<=3 or y>=17″))
# Calculate Bayes factor for each y
# BF = Pr(Y=y | H2) / Pr(Y=y | H1)
y <- 0:20
py.H2 <- 1 / 21 # Exercise: Prove this!
py.H1 <- dbinom(y, 20, .5) # Obviously
BF <- py.H2 / py.H1;
plot(log10(BF) ~ y, type="b")
abline(h=1, col="blue")
abline(h=2, col="red")
legend("top", inset=.10, lty=1, col=c("blue","red"),
legend=c("BF> 10: y<=4 or y>=16″,
“BF>100: y<=2 or y>=18″))