# Illustration of posterior predictive distribution
rm(list=ls()); set.seed(20210128);
# Data are 2 successes in 10 trials
n <- 10; y <- 2; # Prob dist for number of successes in next m trials? m <- 10 # The posterior is theta|y~Beta(a+y, b+n-y) a <- 1; b <- 1 # Uniform prior # Plug in estimate of Pr(y.tilde>0)
1-(1 – y/n)^m
# Approximate the PPD using Monte Carlo
S <- 1e5 # Monte Carlo sample size theta <- rbeta(S, a+y, b+n-y) y.tilde <- rbinom(S, m, theta) mean(y.tilde > 0)
# Compare the dist of y.tilde|theta.hat (“Plug-in”, which
# is the wrong thing) with true dist of y.tilde|y (the
# posterior predictive dist — Right Thing)
PPD <- table(y.tilde)/S; PPD; y.vals <- 0:10 plot(y.vals, dbinom(y.vals,m,y/n), ylab="Probability", xlab="y", cex.lab=1.25, cex.axis=1.25, type="b") lines(y.vals, PPD, type="b", lty=2, pch=19) legend("topright", c("Plug-in","PPD"), lty=c(1,2), pch=c(1,19), bty="n", cex=1.5)