Jim — Sep 11, 2013, 4:41 PM
# Inference about a proportion - part II
# Corrected after class on 9/11/13
# Illustrate
# - several ways of constructing a Bayesian interval estimate
# - prediction
# - sensitivity of the inference / prediction to the choice of prior
# story about hurricanes in Atlantic
# p = proportion that move into land
# illustrate predictive method of choosing a prior
# suppose for a future sample, predict 18% to move on land
# if we see one move on land, predict 22% to move on land
# match this with a beta(3.5, 16) prior
# now observe 20 hurricanes -- 8 move on land
# (I) want a 90% interval estimate for p
# (P) for a future sample of 10, want to predict that at
# least half will reach land
# suppose instead that I used a noninformative uniform prior on p
# how would inferences and prediction change?
# posterior will be Beta(3.5 + 8, 16 + 12)
a <- 3.5 + 8; b <- 16 + 12
# plot posterior
curve(dbeta(x, a, b), 0, 1,
main="Posterior", xlab="P", ylab="Density")
# the triplot shows there is conflict between the prior and data
library(LearnBayes)
triplot(c(3.5, 16), c(8, 12))
# equal-tails 90% interval estimate
b.interval <- qbeta(c(0.05, 0.95), a, b)
b.interval
[1] 0.1799 0.4145
# I decided to use a hpd interval routine in the TeachingDemos package
library(TeachingDemos)
Attaching package: 'TeachingDemos'
The following object is masked from 'package:LearnBayes':
triplot
b.interval2 <- hpd(qbeta, shape1=a, shape2=b, conf=0.90)
b.interval2
[1] 0.1733 0.4067
# compare length of two intervals
diff(b.interval)
[1] 0.2346
diff(b.interval2)
[1] 0.2334
# you should see that the HPD interval is slightly shorter
# prediction -- use LearnBayes function pbetap
pred.prob <- sum(pbetap(c(a, b), 10, 5:10))
pred.prob
[1] 0.1602
# now compare answers using two priors
# with a uniform prior, posterior is Beta(1 + 8, 1 + 12)
a1 <- 1 + 8; b1 <- 1 + 12
# compare "equal-tails" 90% interval estimates
qbeta(c(0.05, 0.95), a, b)
[1] 0.1799 0.4145
qbeta(c(0.05, 0.95), a1, b1)
[1] 0.2450 0.5828
# compare posterior prediction of 5+ hurricanes on land out of 10
sum(pbetap(c(a, b), 10, 5:10))
[1] 0.1602
sum(pbetap(c(a1, b1), 10, 5:10))
[1] 0.4033
# here the inference and the prediction is sensitive to the choice of prior