library(ggplot2)
library(patchwork)
set.seed(536)
theme_set(theme_bw(base_size = 14))
pal <- c(red = "#D62828", blue = "#2E86AB", green = "#06A77D", orange = "#F77F00", dark = "#023047", grey = "#6a6a6a")
find_pdflatex <- function() {
cand <- c(Sys.which("pdflatex"),
path.expand("~/Library/TinyTeX/bin/universal-darwin/pdflatex"),
file.path(Sys.getenv("APPDATA"), "TinyTeX/bin/windows/pdflatex.exe"),
path.expand("~/.TinyTeX/bin/x86_64-linux/pdflatex"))
cand <- cand[nzchar(cand) & file.exists(cand)]
if (!length(cand)) stop("pdflatex not found; install TinyTeX with tinytex::install_tinytex()")
cand[[1]]
}
tikz_png <- function(p, png_file, width = 9, height = 5, dpi = 150) {
for (pkg in c("tikzDevice", "pdftools", "png"))
if (!requireNamespace(pkg, quietly = TRUE)) install.packages(pkg)
options(tikzLatex = find_pdflatex(),
tikzDocumentDeclaration = r"(\documentclass[10pt]{article})",
tikzLatexPackages = c(r"(\usepackage{tikz})",
r"(\usepackage[active,tightpage,psfixbb]{preview})",
r"(\PreviewEnvironment{pgfpicture})",
r"(\setlength\PreviewBorder{0pt})",
r"(\usepackage[T1]{fontenc})", r"(\usepackage{lmodern})",
r"(\usepackage{amsmath})"))
work <- file.path(tempdir(), "tikz_png"); dir.create(work, showWarnings = FALSE)
tex <- file.path(work, "fig.tex"); pdf <- file.path(work, "fig.pdf")
tikzDevice::tikz(tex, width = width, height = height, standAlone = TRUE)
print(p); dev.off()
ok <- system2(getOption("tikzLatex"),
c("-interaction=batchmode", "-halt-on-error",
paste0("-output-directory=", shQuote(work)), shQuote(tex)),
stdout = FALSE, stderr = FALSE)
if (ok != 0 || !file.exists(pdf)) stop("pdflatex failed; see ", file.path(work, "fig.log"))
dir.create(dirname(png_file), recursive = TRUE, showWarnings = FALSE)
png::writePNG(pdftools::pdf_render_page(pdf, page = 1, dpi = dpi), png_file)
png_file
}
# tick labels in math mode, so minus signs and decimals are set by LaTeX
tex_num <- function(v) ifelse(is.na(v), NA,
paste0("$", format(v, scientific = FALSE, drop0trailing = TRUE, trim = TRUE), "$"))
tex_theme <- theme(text = element_text(size = 17), plot.title = element_text(hjust = 0.5),
legend.position = "top", legend.key.width = unit(1.2, "cm"))Lab 7: Bayes’ Rule and Conjugate Priors
STAT 536 — Statistical Computing
Learning Objectives
By the end of today’s lab, you will:
- Compute a posterior with a Bayes table, for two possible values of the unknown
- Derive the Beta–Binomial and Gamma–Poisson posteriors by hand
- Summarize a posterior by its mean, standard deviation, credible interval, and a probability, exactly and by Monte Carlo
- Verify that a posterior mean is a weighted average of the prior mean and the data estimate, and see what more data does to the weights
Setup
Figure labels are typeset by LaTeX. tikz_png() draws a ggplot through tikzDevice, compiles it with pdflatex, and returns a PNG. It needs a LaTeX installation (TinyTeX is enough) and installs tikzDevice, pdftools, and png on first use.
Bayes’ Rule
Recall Bayes’ rule from lecture? The notation is the lecture’s.
| Symbol | Meaning |
|---|---|
| \(\theta\) | the unknown parameter (\(\lambda\) for a Poisson rate) |
| \(\mathbf{x} = (x_1, \ldots, x_n)\) | the observed data |
| \(L(\theta)\) | the likelihood: the model at the observed data, as a function of \(\theta\) |
| \(\pi(\theta)\); \(\pi(\theta \mid \mathbf{x})\) | the prior density (or probability) of \(\theta\); its posterior given the data |
| \(m(\mathbf{x})\) | the normalizing constant: \(\sum_j L(\theta_j)\,\pi(\theta_j)\) for a discrete unknown, \(\int L(\theta)\,\pi(\theta)\,d\theta\) for a continuous one |
| \(\propto\) | equal up to a factor free of \(\theta\) |
\(\pi\) is always a distribution of the unknown. What stands behind the bar says what has been observed.
NoteBayes’ rule
\[\pi(\theta \mid \mathbf{x}) = \frac{L(\theta)\,\pi(\theta)}{m(\mathbf{x})}, \qquad \text{so} \qquad \pi(\theta \mid \mathbf{x}) \propto L(\theta)\,\pi(\theta).\]
Essentially, this says that the posterior is the likelihood times the prior, rescaled to add up (or integrate) to 1.
NoteAlgorithm 1 (Bayes table)
Input: the possible values \(\theta_1, \ldots, \theta_J\), their prior probabilities \(\pi(\theta_j)\), and the data.
- compute the likelihood \(L(\theta_j)\) for each \(j\)
- multiply: \(L(\theta_j)\,\pi(\theta_j)\)
- add the products: \(m(\mathbf{x}) = \sum_j L(\theta_j)\,\pi(\theta_j)\)
- divide each product by \(m(\mathbf{x})\)
Output: the posterior probabilities \(\pi(\theta_j \mid \mathbf{x})\).
Part 1: The Bayes Table
The two coins of lecture. A bag holds a fair coin, \(\theta = 0.5\), and a biased coin, \(\theta = 0.8\). One coin is drawn and flipped \(n = 10\) times. It shows \(k = 7\) heads. The number of heads is \(\Bin(10, \theta)\), so the likelihood is \(L(\theta) = \binom{10}{7}\theta^7(1-\theta)^3\). In R this is dbinom(7, 10, theta).
Exercise 1.1: Two Coins
Implement Algorithm 1.
- Fill in the likelihood and the posterior in
bayes_table(). - Run it with prior probabilities \((0.5, 0.5)\). The coin is drawn at random.
- Run it with prior probabilities \((0.9, 0.1)\). The biased coin is rare.
- The last lines repeat both tables for every possible number of heads, 0 to 10, and keep the posterior probability of the biased coin. Run them as they are.
bayes_table <- function(theta, prior, k, n) {
# TODO: step 1, the likelihood L(theta_j) for each value
likelihood <- ___________
product <- prior * likelihood # step 2
# TODO: steps 3 and 4, the posterior probabilities
posterior <- ___________
data.frame(theta, prior, likelihood, product, posterior)
}
t_equal <- bayes_table(c(0.5, 0.8), c(0.5, 0.5), k = 7, n = 10)
t_rare <- bayes_table(c(0.5, 0.8), c(0.9, 0.1), k = 7, n = 10)
print(t_equal, digits = 3)
print(t_rare, digits = 3)
# the posterior probability of the biased coin, for every number of heads
k_all <- 0:10
biased <- function(prior) {
sapply(k_all, function(k) {
bayes_table(c(0.5, 0.8), prior, k, 10)$posterior[2]
})
}
post_equal <- biased(c(0.5, 0.5))
post_rare <- biased(c(0.9, 0.1))
by_heads <- round(rbind(equal = post_equal, rare = post_rare), 3)
colnames(by_heads) <- k_all
print(by_heads)Report.
- The two tables.
- The posterior probability of the biased coin under each prior.
- The smallest number of heads in 10 flips that makes the biased coin the more probable one, under each prior.
Conjugate Updating
For a continuous unknown the sum in \(m(\mathbf{x})\) is an integral. With a conjugate prior the integral is never computed. The posterior is recognized from its kernel.
NoteKernels
The kernel of a density is the part that depends on the variable.
| Distribution of \(\theta\) | Kernel | Range | Mean | Variance |
|---|---|---|---|---|
| \(\Betad(a, b)\) | \(\theta^{a-1}(1-\theta)^{b-1}\) | \(0 < \theta < 1\) | \(\dfrac{a}{a+b}\) | \(\dfrac{ab}{(a+b)^2(a+b+1)}\) |
| \(\Gam(a, b)\) | \(\theta^{a-1}e^{-b\theta}\) | \(\theta > 0\) | \(\dfrac{a}{b}\) | \(\dfrac{a}{b^2}\) |
\(\Gam(a, b)\) has shape \(a\) and rate \(b\), as in dgamma(x, shape, rate).
NoteAlgorithm 2 (A posterior by hand)
- Write the likelihood \(L(\theta)\). Drop every factor that is free of \(\theta\).
- Multiply by the kernel of the prior.
- Collect the powers and the exponents of \(\theta\). Recognize the kernel and name the posterior, with its parameters.
- Read the summaries off the named distribution: mean, standard deviation, quantiles (
qbeta,qgamma).
NoteDefinition (Credible interval)
A 95% credible interval for \(\theta\) is an interval \([l, u]\) with \[P(l \le \theta \le u \mid \mathbf{x}) = 0.95.\] In this lab \(l\) and \(u\) are the 2.5% and 97.5% quantiles of the posterior.
NoteResult (Weighted average)
For a conjugate pair of this week, \[\text{posterior mean} = w \times \text{prior mean} + (1 - w) \times \text{data estimate}, \qquad w = \frac{n_0}{n_0 + n},\] where \(n_0\) is the number of observations the prior is worth and \(n\) the number of observations in the data. For a \(\Betad(\alpha, \beta)\) prior with binomial data, \(n_0 = \alpha + \beta\) and the data estimate is \(k/n\).
Part 2: Beta–Binomial
The data are \(X \mid \theta \sim \Bin(n, \theta)\), observed \(X = k\). The prior is \(\theta \sim \Betad(\alpha, \beta)\). We use the coin again, \(k = 7\) and \(n = 10\), with two priors.
- The flat prior \(\Betad(1, 1)\). Every \(\theta\) between 0 and 1 is equally plausible before the data.
- The fair-coin prior \(\Betad(5, 5)\). It is centered at \(0.5\).
Exercise 2.1: The Posterior by Hand
Follow Algorithm 2 on paper.
- Write the likelihood \(L(\theta)\). Drop every factor that is free of \(\theta\).
- Multiply by the kernel of the \(\Betad(\alpha, \beta)\) prior.
- Collect the powers of \(\theta\) and of \(1 - \theta\). Name the posterior and give its two parameters.
Then fill in beta_update() with your two parameters and run the code chunk.
Report.
- Your derivation, in three lines.
- The posterior parameters under each prior.
Exercise 2.2: Summaries and a Picture
- Fill in
beta_summary(): the two ends of the 95% credible interval, and \(P(\theta > 0.5 \mid \mathbf{x})\). The mean and the standard deviation are given. - Run it for both posteriors.
- Run the figure chunk. It draws the prior, the likelihood rescaled to integrate to 1, and the posterior, for each prior. The dotted line is the observed proportion \(k/n = 0.7\).
beta_summary <- function(a, b) {
# TODO: the 2.5% and 97.5% quantiles, and P(theta > 0.5 | x)
c(mean = a / (a + b),
sd = sqrt(a * b / ((a + b)^2 * (a + b + 1))),
lower = ___________, upper = ___________,
p_gt_half = ___________)
}
sum_flat <- beta_summary(post_flat[["alpha"]], post_flat[["beta"]])
sum_fair <- beta_summary(post_fair[["alpha"]], post_fair[["beta"]])
print(round(rbind(flat = sum_flat, fair = sum_fair), 3))three_curves <- function(grid, prior, lik, post, xlab, title) {
lev <- c("prior", "likelihood (rescaled)", "posterior")
curve <- rep(factor(lev, levels = lev), each = length(grid))
df <- data.frame(x = rep(grid, 3), density = c(prior, lik, post), curve)
ggplot(df, aes(x, density, color = curve, linetype = curve,
linewidth = curve, alpha = curve)) +
geom_line() +
scale_color_manual(values = unname(pal[c("blue", "orange", "red")]),
name = NULL) +
scale_linetype_manual(values = c("dashed", "solid", "solid"),
name = NULL) +
scale_linewidth_manual(values = c(0.8, 2.4, 1), name = NULL) +
scale_alpha_manual(values = c(1, 0.55, 1), name = NULL) +
scale_x_continuous(labels = tex_num) +
scale_y_continuous(labels = tex_num) +
labs(x = xlab, y = "density", title = title) + tex_theme
}
theta <- seq(0, 1, length.out = 400)
# the likelihood theta^7 (1 - theta)^3, rescaled to integrate to 1,
# is the Beta(8, 4) density
coin_panel <- function(a0, b0, title) {
pn <- beta_update(a0, b0, k = 7, n = 10)
three_curves(theta, dbeta(theta, a0, b0), dbeta(theta, 8, 4),
dbeta(theta, pn[["alpha"]], pn[["beta"]]),
r"($\theta$)", title) +
geom_vline(xintercept = 0.7, color = pal[["dark"]],
linetype = "dotted") +
coord_cartesian(ylim = c(0, 3.8))
}
p <- coin_panel(1, 1, r"(Prior $\mathrm{Beta}(1, 1)$)") +
coin_panel(5, 5, r"(Prior $\mathrm{Beta}(5, 5)$)") +
plot_layout(guides = "collect") & theme(legend.position = "top")
knitr::include_graphics(tikz_png(p, "figures/lab07/coin-posterior.png",
width = 10, height = 4.8))Report.
- The table of summaries.
- The figure.
- Which prior gives the narrower credible interval, and why, in one sentence.
Exercise 2.3: The Weighted Average, and Two Batches
- Fill in the weight \(w\) and the weighted average in
weighted_mean(). Use the Result box: for a \(\Betad(\alpha, \beta)\) prior, \(n_0 = \alpha + \beta\). - Compare the weighted average with the posterior mean of Exercise 2.2, for both priors.
- Update in two batches. Start from the flat prior’s posterior after the first 10 flips. Update it on 30 more flips with 21 heads. Fill in the call that does this.
- Compare with one update of the flat prior on all 40 flips, 28 heads.
weighted_mean <- function(alpha, beta, k, n) {
# TODO: the weight on the prior mean
w <- ___________
# TODO: w * prior mean + (1 - w) * data estimate
c(w = w, weighted = ___________)
}
wm_flat <- weighted_mean(1, 1, 7, 10)
wm_fair <- weighted_mean(5, 5, 7, 10)
print(round(rbind(flat = c(wm_flat, direct = sum_flat[["mean"]]),
fair = c(wm_fair, direct = sum_fair[["mean"]])), 4))
# TODO: update post_flat on 21 heads in 30 more flips
two_steps <- ___________
one_step <- beta_update(1, 1, k = 28, n = 40)
print(rbind(two_steps, one_step))Report.
- The weight \(w\) and both versions of the posterior mean, for each prior.
- The two-batch posterior and the one-batch posterior.
- How many flips each prior is worth.
Exercise 2.4: Summaries by Monte Carlo
Recall Monte Carlo integration from week 1? A posterior with a name can be sampled, and every summary is then an average over the draws. The posterior under the flat prior is \(\Betad(8, 4)\). Draw \(S = 10^5\) values from it.
- Estimate \(P(\theta > 0.5 \mid \mathbf{x})\) by the fraction of draws above \(0.5\). Give its Monte Carlo standard error.
- Compute the odds \(\theta/(1 - \theta)\) for each draw.
- The code chunk then gives the posterior median and the 95% credible interval of the odds, from the draws and exactly. The odds are increasing in \(\theta\), so a quantile of the odds is the same quantile of \(\theta\) pushed through \(q/(1 - q)\).
set.seed(536)
S <- 1e5
draws <- rbeta(S, 8, 4)
# TODO: the fraction of draws above 0.5, and its Monte Carlo standard error
p_mc <- ___________
p_se <- ___________
# TODO: the odds, draw by draw
odds <- ___________
odds_mc <- quantile(odds, c(0.025, 0.5, 0.975))
q <- qbeta(c(0.025, 0.5, 0.975), 8, 4)
odds_exact <- q / (1 - q) # the quantiles of theta, transformed
print(round(c(p_mc = p_mc, p_se = p_se, p_exact = 1 - pbeta(0.5, 8, 4)), 4))
print(round(rbind(mc = odds_mc, exact = odds_exact), 3))Report.
- The Monte Carlo estimate of \(P(\theta > 0.5 \mid \mathbf{x})\), its standard error, and the exact value.
- The median and the 95% credible interval of the odds, from the draws and exactly.
- How many standard errors the estimate is from the exact value.
Part 3: Gamma–Poisson
Recall the Poisson counts of the accept-reject lab? Ten counts with an unknown rate \(\lambda\) and a \(\Gam(2, 1)\) prior. That lab stated the posterior and sampled from it. Here we derive it.
The data are \(X_1, \ldots, X_n \mid \lambda \iid \Pois(\lambda)\), so each count has probability \(\lambda^{x_i}e^{-\lambda}/x_i!\). The prior is \(\lambda \sim \Gam(\alpha, \beta)\), with density \(\dfrac{\beta^\alpha}{\Gamma(\alpha)}\lambda^{\alpha-1}e^{-\beta\lambda}\) and mean \(\alpha/\beta\).
Exercise 3.1: The Posterior and Its Mean
On paper:
- Write the likelihood \(L(\lambda)\) of the \(n\) counts. Drop every factor that is free of \(\lambda\).
- Multiply by the kernel of the prior.
- Name the posterior and give its two parameters.
- Write the posterior mean as a weighted average of the prior mean \(\alpha/\beta\) and the sample mean \(\bar x\). Identify the weight \(w\) on the prior mean.
Then fill in gamma_update() and the weight w, and run the code chunk and the figure chunk. The dotted line in the figure is \(\bar x\).
x <- c(3, 5, 4, 6, 3, 4, 5, 7, 4, 3) # the counts of the accept-reject lab
gamma_update <- function(alpha, beta, x) {
# TODO: the posterior shape and rate from your derivation
c(alpha = ___________, beta = ___________)
}
post <- gamma_update(2, 1, x)
a_n <- post[["alpha"]]; b_n <- post[["beta"]]
gp_summary <- c(mean = a_n / b_n, sd = sqrt(a_n) / b_n,
lower = qgamma(0.025, a_n, b_n),
upper = qgamma(0.975, a_n, b_n))
# TODO: the weight on the prior mean, for the Gamma(2, 1) prior
w <- ___________
prior_mean <- 2 / 1 # alpha / beta of the Gamma(2, 1) prior
weighted <- w * prior_mean + (1 - w) * mean(x)
print(post)
print(round(gp_summary, 3))
print(round(c(w = w, weighted = weighted, direct = a_n / b_n), 4))lambda <- seq(0, 8, length.out = 400)
# the likelihood lambda^sum(x) exp(-n lambda), rescaled to integrate to 1,
# is the Gamma(sum(x) + 1, n) density
p <- three_curves(lambda, dgamma(lambda, 2, 1),
dgamma(lambda, sum(x) + 1, length(x)),
dgamma(lambda, a_n, b_n), r"($\lambda$)",
"Prior, likelihood, and posterior for the Poisson rate") +
geom_vline(xintercept = mean(x), color = pal[["dark"]],
linetype = "dotted")
knitr::include_graphics(tikz_png(p, "figures/lab07/gamma-poisson.png",
width = 9, height = 4.5))Report.
- Your derivation: the four steps.
- The posterior parameters, mean, standard deviation, and 95% credible interval.
- The weight on the prior mean, and both versions of the posterior mean.
- The figure.
Exercise 3.2: More Data, Two Priors
The code chunk below is complete. Run it as is. It simulates 1000 counts from \(\Pois(4)\), so the true rate is \(4\). It updates two priors on the first 10 counts, the first 100, and all 1000.
- \(\Gam(2, 1)\): prior mean \(2\), worth 1 observation.
- \(\Gam(50, 5)\): prior mean \(10\), worth 5 observations.
Then run the figure chunk. It draws the two posteriors at each \(n\). The dotted line is the true rate.
set.seed(536)
y <- rpois(1000, 4) # the true rate is 4
one_row <- function(m, prior) {
p <- gamma_update(prior[1], prior[2], y[1:m])
data.frame(n = m, alpha0 = prior[1], beta0 = prior[2],
xbar = mean(y[1:m]), w = prior[2] / (prior[2] + m),
mean = p[["alpha"]] / p[["beta"]],
lower = qgamma(0.025, p[["alpha"]], p[["beta"]]),
upper = qgamma(0.975, p[["alpha"]], p[["beta"]]))
}
grow <- do.call(rbind, lapply(c(10, 100, 1000), function(m) {
rbind(one_row(m, c(2, 1)), one_row(m, c(50, 5)))
}))
print(round(grow, 3))lambda <- seq(2, 8, length.out = 500)
dens <- do.call(rbind, lapply(seq_len(nrow(grow)), function(i) {
a0 <- grow$alpha0[i]; b0 <- grow$beta0[i]; m <- grow$n[i]
p <- gamma_update(a0, b0, y[1:m])
data.frame(lambda, n = m,
density = dgamma(lambda, p[["alpha"]], p[["beta"]]),
prior = sprintf("prior $\\mathrm{Gamma}(%d, %d)$", a0, b0))
}))
panels <- sprintf("$n = %d$", c(10, 100, 1000))
dens$panel <- factor(sprintf("$n = %d$", dens$n), levels = panels)
p <- ggplot(dens, aes(lambda, density, linetype = prior)) +
geom_line(color = pal[["red"]], linewidth = 1) +
geom_vline(xintercept = 4, color = pal[["dark"]], linetype = "dotted") +
facet_wrap(~ panel, scales = "free_y") +
scale_linetype_manual(values = c("dashed", "solid"), name = NULL) +
scale_x_continuous(labels = tex_num) +
scale_y_continuous(labels = tex_num) +
labs(x = r"($\lambda$)", y = "density",
title = r"(Posteriors under two priors as $n$ grows)") +
tex_theme
knitr::include_graphics(tikz_png(p, "figures/lab07/more-data.png",
width = 10, height = 4.2))Report.
- The table.
- The figure.
- The distance between the two posterior means at each \(n\).
- Whether each 95% credible interval contains the true rate \(4\), at each \(n\).
Summary and Reflection
Takeaways
- One table, any prior. Bayes’ rule is a multiplication and a division: prior times likelihood, over the sum.
- Conjugate updating is addition. The posterior parameters are the prior parameters plus counts from the data. Two batches and one batch give the same posterior.
- The posterior mean is a weighted average. The weight on the prior is its sample size over the total.
- Quantiles for a named posterior, draws for any posterior.
qbetaandqgammagive credible intervals. Draws give every other summary, each with a Monte Carlo standard error. - Data swamp the prior. How fast depends on how many observations the prior is worth.
Self-Assessment Questions
- In the two-coin table, which column changes when the prior changes, and which does not?
- Why can the factor \(\binom{n}{k}\) be dropped from the likelihood before multiplying by the prior?
- A \(\Betad(2, 8)\) prior is updated with 3 successes in 10 trials. Give the posterior, its mean, and the weight on the prior mean, without a computer.
- For Poisson data with a \(\Gam(\alpha, \beta)\) prior, what are the two ingredients of the posterior mean, and what happens to the weight on the prior as \(n\) grows?
- A colleague reports a 95% credible interval and calls it a 95% confidence interval. State what each 95% means, in one sentence each.
- You have \(10^5\) draws from a posterior. How do you estimate the posterior probability that the parameter exceeds 3, and how accurate is the estimate?