# Chihara and Hesterberg, chap 10 # Chapter 10 Bayesian Methods # R scripts #------------------------- # Section 10.2 binomial data, discrete prior distributions #----------------------------------- # Example 10.1 tennis theta <- seq(0, 1, by = .1) prior <- c(0, .02, .03, .05, .1, .15, .2, .25, .15, .05, 0) # results : WLL likelihood <- theta * (1 - theta)^2 constant <- sum(prior * likelihood) posterior <- prior * likelihood / constant mean.prior <- sum(theta * prior); mean.prior # prior mean mean.posterior1 <- sum(theta * posterior); mean.posterior1 # posterior mean # more results : W 2, L 3 # use total results : W 3, L5 likelihood2 <- theta^3 * (1 - theta)^5 # 3 success, 5 fail constant2 <- sum(prior * likelihood2) posterior2 <- prior * likelihood2 / constant2 # use additional results : W 2, L 3 likelihood3 <- theta^2 * (1 - theta)^3 constant3 <- sum(posterior * likelihood3) posterior3 <- posterior * likelihood3 / constant3 posterior3 # not shown, matches posterior2 mean.posterior2 <- sum(theta*posterior2); mean.posterior2 # posterior mean plot(theta, prior, type = "b", col="palegreen4", ylim = c(0, max(posterior3)), ylab = "probability", main="Tennis") lines(theta, posterior, type = "b", lty = 2, col="paleturquoise3") lines(theta, posterior2, type = "b", lty = 3, col="palevioletred3") legend("topleft", legend = c("prior", "posterior1", "posterior2"), col=c("palegreen4", "paleturquoise3", "palevioletred3"), lty = 1:3) data.frame(mean.prior, mean.posterior1, mean.posterior2) #----------------------------------- # Example 10.2 exit poll theta <- 4:6 / 10 prior <- c(.35, .50, .15) # results : CCM likelihood <- theta^2 * (1 - theta) constant <- sum(prior * likelihood) posterior <- prior * likelihood / constant; posterior # more results : CMCC likelihood2 <- theta^3 * (1 - theta) constant2 <- sum(posterior * likelihood2) posterior2 <- posterior * likelihood2 / constant2; posterior2 data.frame(theta, prior, posterior, posterior2) plot(theta, prior, type = "b", col="palegreen4", ylim = c(0, max(posterior2)), ylab = "probability", main="Exit Poll") lines(theta, posterior, type = "b", lty = 2, col="paleturquoise3") lines(theta, posterior2, type = "b", lty = 3, col="palevioletred3") legend("topleft", legend = c("prior", "posterior", "posterior2"), col=c("palegreen4", "paleturquoise3", "palevioletred3"), lty = 1:3) #------------------------- # Section 10.3 binomial data, continuous prior distributions #------------------------- # beta distributions library(RColorBrewer) display.brewer.all() palette.GnBu <- brewer.pal(9, "GnBu") curve(dbeta(x, 0.5, 6), from=0, to=1, col=palette.GnBu[[4]], ylim=c(0,8), xlab=expression(theta), ylab = expression(paste("density of Beta(", alpha, ", ", beta,")")), main=expression(paste("Beta(", alpha, ", ", beta,")"))) curve(dbeta(x, 1.8, 6), col=palette.GnBu[[5]], add=TRUE) curve(dbeta(x, 1, 1), col=palette.GnBu[[6]], add=TRUE) curve(dbeta(x, 6, 6), col=palette.GnBu[[7]], add=TRUE) curve(dbeta(x, 30, 30), col=palette.GnBu[[8]], add=TRUE) curve(dbeta(x, 5, 1), col=palette.GnBu[[9]], add=TRUE) legend("topright", legend = c("beta(0.5, 6)", "beta(1.8, 6)", "beta(1, 1)", "beta(6, 6)", "beta(30, 30)", "beta(5, 1)"), lwd=2, col=palette.GnBu[4:9]) #------------------------- # Example 10.3 Jon Stewert n <- 200 x <- 110 # frequentist theta.hat <- x / n; thetaHat # Agresti-Coull confidence interval (p.183) alpha <- .05 n.tilda <- n + 4 theta.tilda <- (x + 2) / n.tilda z.a <- qnorm(1 - alpha/2) e <- z.a * sqrt(theta.tilda * (1 - theta.tilda) / n.tilda) frequentist.ci <- theta.tilda + c(-e, e); frequentist.ci # => 0.4807377 0.6173015 # Bayesian A # prior density : pi(theta) == beta(1, 1) # posterior density : p(theta | X === 110) == beta(1 + 110, 1 + 200 - 110) == beta(111, 91) display.brewer.all() palette.RdPu <- brewer.pal(9, "RdPu") alpha <- 111 beta <- 91 curve(dbeta(x, alpha, beta), from=0, to=1, lty=1, col=palette.RdPu[[7]], ylim=c(0,11.5), xlab=expression(x), ylab = expression(paste("density of Beta(", alpha, ", ", beta,")")), main="Jon Stewert : A") curve(dbeta(x, 1, 1), lty=2, col=palette.RdPu[[5]], add=TRUE) legend(.05, 11, legend = c("prior", "posterior"), lwd=2, lty=1:2, col=c(palette.RdPu[c(5, 7)])) posterior.mean <- alpha / (alpha + beta); posterior.mean # credible interval bayesian.A.ci <- c(qbeta(.025, alpha, beta), qbeta(.975, alpha, beta)); bayesian.A.ci # => 0.4806705 0.6174106 # P(theta >= 0.5 | X == 110) # => 0.920917 (Mathematica) 1 - pbeta(.5, alpha, beta) # => 0.9209173 # Bayesian B n <- 200 x <- 110 # prior mean <- .58 sd <- .03 # from Mathematica alpha.B <- 156.407 beta.B <- 113.26 # posterior a <- x + alpha.B b <- n - x + beta.B # beta(a, b) == beta(110 + alpha.B, 200 - 110 + beta.B) == beta(266.407, 203.26) posterior.mean <- a / (a + b); posterior.mean # => 0.5672253 posterior.var <- a * b / ((a + b)^2 * (a + b + 1)); posterior.var posterior.sd <- sqrt(posterior.var); posterior.sd # => 0.02283767 # credible interval bayesian.B.ci <- c(qbeta(.025, a, b), qbeta(.975, a, b)); bayesian.B.ci # => 0.5222139 0.6116944 alpha <- a beta <- b curve(dbeta(x, alpha, beta), from=0, to=1, lty=1, col=palette.RdPu[[8]], ylim=c(0,18), xlab=expression(x), ylab = expression(paste("density of Beta(", alpha, ", ", beta,")")), main="Jon Stewert : B") curve(dbeta(x, alpha.B, beta.B), lty=2, col=palette.RdPu[[6]], add=TRUE) legend(.05, 17, legend = c("prior", "posterior"), lwd=2, lty=1:2, col=c(palette.RdPu[6], palette.RdPu[8])) #------------------------- # Example 10.4 Continuous data : trout # fish length X ~ N(mu, 8^2) # prior distribution mu ~ N(50, 6^2) n <- 15 xBar <- 45 sigma <- 8 mu0 <- 50 sigma0 <- 6 precision.prior <- 1 / sigma0^2 precision.data <- n / sigma^2 mu1 <- (precision.prior * mu0 + precision.data * xBar) / (precision.prior + precision.data); mu1 # => 45.5298 precision.posterior <- precision.prior + precision.data; precision.posterior var.posterior <- 1 / precision.posterior; var.posterior sd.posterior <- sqrt(var.posterior); sd.posterior # => 1.953092 # posterior : mu | data ~ N(mu1, sd.posterior^2) == N(45.5298, 1.953092^2) library(RColorBrewer) display.brewer.all() palette.OrRd <- brewer.pal(9, "OrRd") help("points") # available characters curve(dnorm(x, mu0, sigma0), from=35, to=65, lty=2, col=palette.OrRd[[4]], ylim=c(0, .20), xlab=expression(x), ylab = "density for normal", main="Trout") curve(dnorm(x, mu1, sd.posterior), lty=1, col=palette.OrRd[[8]], add=TRUE) points(mu1, 0, pch=8, cex=1.4, col="tomato4") legend(53, .17, legend = c("prior", "posterior"), lwd=2, lty=c(2, 1), col=palette.OrRd[c(4, 8)]) #------------------------- # Chapter 10.5 Sequential data #------------------------- # Example 10.5 poll of students # assume a flat prior a0 <- 1 b0 <- 1 # first poll x1 <- 10 n1 <- 38 a1 <- a0 + x1; a1 b1 <- b0 + n1 - x1; b1 # posterior : P(theta | x == 10) == beta(a + x, b + n - x) == beta(11, 29) # E(theta | x == 10) == 11 / (11 + 29) == a1 / (a1 + b1) mu1 <- a1 / (a1 + b1); mu1 # == (a0 + x1) / (a0 + b0 + n1) # => 0.275 # second poll a1 <- 11 b1 <- 29 x2 <- 5 n2 <- 15 a2 <- a1 + x2; a2 b2 <- b1 + n2 - x2; b2 # posterior : P(theta | x == 5) == beta(a1 + x2, b1 + n2 - x2) == beta(16, 39) mu2 <- a2 / (a2 + b2); mu2 # => 0.2909091 # E(theta | x2 == 5) == 16 / (16 + 39) # all the data x <- x1 + x2 n <- n1 + n2 a <- a0 + x; a b <- b0 + n - x; b # posterior : P(theta | x == 15) == beta(16, 39) #------------------------- # Example web pages n <- c(1874, 1867, 1871, 1868, 1875, 1875) X <- c(52, 41, 55, 49, 39, 39) alpha <- X # vector of posterior parameters beta <- n - X # vector of posterior parameters N <- 10^5 # replications theta <- matrix(0.0, nrow = N, ncol = 6) for (j in 1:6) { theta[, j] <- rbeta(N, alpha[j], beta[j]) } probBest <- numeric(6) # vector for results best <- apply(theta, 1, max) # maximum of each row (margin=1 is code for rows) for (j in 1:6) { probBest[j] <- mean(theta[, j] == best) } probBest plot(theta[1:10^4, 1], theta[1:10^4, 3], col="orange3", pch = ".", main="Sequential Data : 1 vs 3") abline(0, 1, col="dark red") text(.037, .042, substitute(theta[3] > theta[1])) text(.042, .037, substitute(theta[1] > theta[3])) plot(theta[1:10^4, 2], theta[1:10^4, 3], col="orange3", pch = ".", main="Sequential Data : 2 vs 3") abline(0, 1, col="dark red") text(.034, .040, substitute(theta[3] > theta[2])) text(.034, .030, substitute(theta[2] > theta[3])) #-------------------------