Jim — Sep 6, 2013, 3:20 PM
################################################
# R demonstration September 6, 2013
# Learning About a Proportion
################################################
# load LearnBayes package (assuming it has been installed)
library(LearnBayes)
# using a discrete prior
# suppose p has values .1, .3, .5, .7, .9
# with weights (probs) 4, 10, 4, 2, 1
# set up prior
p <- seq(0.1, 0.9, .2)
p
[1] 0.1 0.3 0.5 0.7 0.9
prior <- c(4, 10, 4, 2, 1)
prior <- prior / sum(prior)
prior
[1] 0.19048 0.47619 0.19048 0.09524 0.04762
# plot prior
plot(p, prior, type="h")
# observe some data: sample of 20 with 4 successes, 16 failures
s <- 4; f <- 16
# compute posterior probabilities
?pdisc
starting httpd help server ... done
pdisc(p, prior, c(s, f)) -> post.prob
# display and plot
post.prob
[1] 2.135e-01 7.755e-01 1.099e-02 5.955e-06 1.890e-13
plot(p, post.prob, type="h",
main="Plot of Posterior Probabilities")
# continuous prior
# suppose I believe P(p < .2) = 0.5, P(p < .4) = 0.9
# want to find a beta(a, b) prior which matches these beliefs
quantile1 <- list(p=.5, x=.2)
quantile2 <- list(p=.9, x=.4)
ab <- beta.select(quantile1, quantile2)
ab
[1] 2.07 7.32
# plot a beta prior
curve(dbeta(x, ab[1], ab[2]), 0, 1,
main="Beta(2.07, 7.32) prior")
# Recall we observed data.
s
[1] 4
f
[1] 16
# illustrate a Bayesian triplot for a proportion
triplot(ab, c(s, f))
# plot the beta posterior
curve(dbeta(x, ab[1] + s,
ab[2] + f), 0, 1,
main="Beta(6.07, 23.32) Posterior")
# illustrate summarizing this posterior
ab1 <- c(ab[1] + s, ab[2] + f)
# P(theta < 0.2)
pbeta(0.2, ab1[1], ab1[2])
[1] 0.501
# P( 0.1 < theta < 0.2)
diff(pbeta(c(0.1, 0.2),
ab1[1], ab1[2]))
[1] 0.4467
# interval that brackets 90% of the probability
round(qbeta(c(0.05, 0.95),
ab1[1], ab1[2]), 3)
[1] 0.098 0.338
# illustrate computing the (prior) predictive distribution
# future sample of size 20, what is the chance of at least 5 successes?
# illustrate learning about p
# using simulation
ab1
[1] 6.07 23.32
rbeta(1000, ab1[1], ab1[2]) ->
sim.p
hist(sim.p)
plot(density(sim.p))
# add real beta posterior
curve(dbeta(x, ab1[1], ab1[2]),
add=TRUE, col="red", lwd=3)
# find a 90% interval estimate
quantile(sim.p, c(0.05, 0.95))
5% 95%
0.1046 0.3443
# what's P(p > 0.6)?
sum(sim.p > .6) / 1000
[1] 0
# exact value is
1 - pbeta(0.6, ab1[1], ab1[2])
[1] 5.335e-06
################################################