Bayesian Statistical Inference for Basic Science: Part 2
From general statistical inference to posterior learning, uncertainty, and scientific data analysis
Dr. Debashis Chatterjee
10 July 2026
1 1. Why this part comes after general statistical inference
In Part 1, we learned the language of general statistical inference: sample, estimator, standard error, confidence interval, hypothesis test, likelihood, regression, bootstrap, and model checking.
This Part 2 keeps the same scientific goal but changes the inferential viewpoint.
Instead of saying only \[ \widehat{\theta}\quad\text{estimates}\quad \theta, \] Bayesian inference treats the unknown scientific quantity itself as uncertain and asks for \[ \pi(\theta\mid y), \] that is, the probability distribution of the unknown parameter after observing the data.
For basic science researchers, the key idea is this:
Bayesian inference is probabilistic inference about unknown scientific quantities, conditional on the observed data and the assumed scientific-statistical model.
The goal is not merely to replace a p-value by a prior. The goal is to express scientific uncertainty, combine data with structure, and produce uncertainty-aware predictions and decisions.
2 2. Probability environment and Bayesian statistical model
2.1 2.1 Data-generating probability space
A probability model begins with a probability space \[ (\Omega,\mathcal F,P). \]
Here:
- \(\Omega\) is the sample space of possible outcomes.
- \(\mathcal F\) is a sigma-field of events.
- \(P\) is a probability measure.
A random variable is a measurable map \[ Y:\Omega\to\mathcal Y. \]
In statistical inference, we usually do not know the true probability law. Instead, we specify a family of possible laws \[ \mathcal P=\{P_\theta:\theta\in\Theta\}, \] where \(\theta\) is the unknown parameter.
Examples:
| Scientific situation | Observation \(Y\) | Parameter \(\theta\) | Statistical model |
|---|---|---|---|
| Photon counting | source counts | source intensity | Poisson |
| Detection experiment | success/failure | detection probability | Binomial |
| Laboratory measurement | noisy reading | true mean signal | Normal |
| River flow series | annual discharge | changepoint, means, variance | time-series/changepoint |
| Galaxy velocities | velocities | mixture components | mixture model |
2.2 2.2 Bayesian layer: prior and joint model
Bayesian inference adds a probability distribution on the parameter space: \[ \theta\sim\Pi. \]
If a density exists, we write \[ \pi(\theta). \]
The data model is \[ Y\mid \theta \sim P_\theta, \] or, in density form, \[ p(y\mid \theta). \]
Together these define the joint model \[ p(y,\theta)=p(y\mid\theta)\pi(\theta). \]
The posterior is the conditional distribution \[ \pi(\theta\mid y)=\frac{p(y\mid\theta)\pi(\theta)}{m(y)}, \] where \[ m(y)=\int_\Theta p(y\mid\theta)\pi(\theta)d\theta \] is the marginal likelihood, also called evidence.
2.3 2.3 What does the posterior mean?
After observing data \(y\), the posterior \(\pi(\theta\mid y)\) tells us which parameter values remain plausible under the model.
From the posterior we can compute:
- posterior mean: \(E(\theta\mid y)\),
- posterior median,
- credible interval,
- posterior tail probability, such as \(P(\theta>0\mid y)\),
- posterior predictive distribution: \[ p(\widetilde y\mid y)=\int p(\widetilde y\mid\theta)\pi(\theta\mid y)d\theta. \]
3 3. A first complete example: Binomial detection experiment
Imagine a telescope detector. During \(n\) independent trials, the detector either detects a weak signal or it does not.
Let \[ Y\mid\theta\sim\text{Binomial}(n,\theta), \] where \(\theta\) is the unknown detection probability.
Use a Beta prior: \[ \theta\sim\text{Beta}(a,b). \]
Then the posterior is \[ \theta\mid Y=y\sim\text{Beta}(a+y,b+n-y). \]
This is conjugacy: prior and posterior belong to the same family.
3.1 3.1 Posterior concentration as larger data arrive
Let the true detection probability be \(\theta_0=0.37\). We will reveal the data sequentially and watch the posterior concentrate.
theta_true <- 0.37
nmax <- 1200
yseq <- rbinom(nmax, size = 1, prob = theta_true)
n_grid <- c(5, 10, 20, 50, 100, 250, 500, 1200)
a0 <- 2; b0 <- 2
xgrid <- seq(0.001, 0.999, length.out = 900)
dens_list <- list()
sum_tab <- data.frame()
for(n in n_grid) {
y <- sum(yseq[1:n])
aa <- a0 + y
bb <- b0 + n - y
dens_list[[as.character(n)]] <- data.frame(n = factor(n, levels = n_grid),
theta = xgrid,
density = dbeta(xgrid, aa, bb))
ci <- qbeta(c(0.025, 0.975), aa, bb)
sum_tab <- rbind(sum_tab, data.frame(n = n,
y = y,
sample_prop = y/n,
post_mean = aa/(aa+bb),
ci_low = ci[1],
ci_high = ci[2],
width = ci[2]-ci[1]))
}
post_dens <- do.call(rbind, dens_list)
pretty_table(sum_tab, digits = 4,
caption = "Sequential Bayesian learning for a Binomial detection probability.")| n | y | sample_prop | post_mean | ci_low | ci_high | width |
|---|---|---|---|---|---|---|
| 5 | 0 | 0.0000 | 0.2222 | 0.0319 | 0.5265 | 0.4947 |
| 10 | 1 | 0.1000 | 0.2143 | 0.0504 | 0.4545 | 0.4041 |
| 20 | 4 | 0.2000 | 0.2500 | 0.1023 | 0.4370 | 0.3347 |
| 50 | 16 | 0.3200 | 0.3333 | 0.2152 | 0.4632 | 0.2480 |
| 100 | 37 | 0.3700 | 0.3750 | 0.2849 | 0.4697 | 0.1848 |
| 250 | 84 | 0.3360 | 0.3386 | 0.2818 | 0.3978 | 0.1161 |
| 500 | 173 | 0.3460 | 0.3472 | 0.3063 | 0.3893 | 0.0830 |
| 1200 | 445 | 0.3708 | 0.3713 | 0.3442 | 0.3987 | 0.0545 |
ggplot(post_dens, aes(theta, density, color = n)) +
geom_line(linewidth = 1) +
geom_vline(xintercept = theta_true, linetype = 2, color = "black") +
labs(title = "Posterior concentration as more detector trials arrive",
subtitle = "The dashed vertical line is the true detection probability used in simulation.",
x = expression(theta), y = "Posterior density", color = "n")
3.2 3.2 Credible interval shrinkage
ggplot(sum_tab, aes(n, post_mean)) +
geom_ribbon(aes(ymin = ci_low, ymax = ci_high), alpha = 0.18) +
geom_line(linewidth = 1.1) +
geom_point(size = 2.5) +
geom_hline(yintercept = theta_true, linetype = 2) +
scale_x_log10(breaks = n_grid) +
labs(title = "Posterior mean and 95% credible interval as information grows",
subtitle = "The interval shrinks, and the posterior mean approaches the true value in this correctly specified simulation.",
x = "Sample size n, log scale", y = expression(theta))
ggplot(sum_tab, aes(n, width)) +
geom_line(linewidth = 1.1, color = "firebrick") +
geom_point(size = 2.5, color = "firebrick") +
scale_x_log10(breaks = n_grid) +
labs(title = "Posterior uncertainty decreases with increasing information",
x = "Sample size n, log scale", y = "Width of 95% credible interval")
4 4. Simulation verification: posterior coverage and learning
A single simulated dataset is illustrative. A researcher should also ask: if the same experiment were repeated many times, how often do Bayesian credible intervals cover the true value?
set.seed(101)
R <- 500
nvals <- c(20, 50, 100, 250, 500)
cover_tab <- data.frame()
for(n in nvals) {
covered <- numeric(R)
width <- numeric(R)
err <- numeric(R)
for(r in seq_len(R)) {
y <- rbinom(1, n, theta_true)
aa <- a0 + y
bb <- b0 + n - y
ci <- qbeta(c(0.025, 0.975), aa, bb)
pm <- aa/(aa + bb)
covered[r] <- as.numeric(theta_true >= ci[1] && theta_true <= ci[2])
width[r] <- ci[2] - ci[1]
err[r] <- abs(pm - theta_true)
}
cover_tab <- rbind(cover_tab, data.frame(n = n,
coverage = mean(covered),
mean_width = mean(width),
mean_abs_error = mean(err)))
}
pretty_table(cover_tab, digits = 4,
caption = "Repeated-experiment verification under a correctly specified Beta-Binomial model.")| n | coverage | mean_width | mean_abs_error |
|---|---|---|---|
| 20 | 0.968 | 0.3721 | 0.0732 |
| 50 | 0.966 | 0.2523 | 0.0509 |
| 100 | 0.958 | 0.1839 | 0.0366 |
| 250 | 0.944 | 0.1184 | 0.0251 |
| 500 | 0.932 | 0.0842 | 0.0185 |
cover_long <- rbind(
data.frame(n = cover_tab$n, quantity = "Coverage", value = cover_tab$coverage),
data.frame(n = cover_tab$n, quantity = "Mean interval width", value = cover_tab$mean_width),
data.frame(n = cover_tab$n, quantity = "Mean absolute error", value = cover_tab$mean_abs_error)
)
ggplot(cover_long, aes(n, value, color = quantity)) +
geom_line(linewidth = 1.1) + geom_point(size = 2.5) +
scale_x_log10(breaks = nvals) +
facet_wrap(~ quantity, scales = "free_y") +
labs(title = "Bayesian learning verified by repeated simulation",
subtitle = "Correct model: posterior error and interval width decrease as n increases.",
x = "Sample size n, log scale", y = "Value", color = "") +
theme(legend.position = "none")
5 5. Prior influence and why data eventually dominate
Bayesian inference does not say that the prior is irrelevant. It says prior and likelihood are combined. With increasing data, a regular prior usually becomes less influential.
Consider \[ Y_i\mid\mu\sim N(\mu,\sigma^2),\qquad \sigma^2\text{ known}, \] and prior \[ \mu\sim N(m_0,s_0^2). \]
Then \[ \mu\mid y_{1:n}\sim N(m_n,s_n^2), \] where \[ s_n^2=\left(\frac{1}{s_0^2}+\frac{n}{\sigma^2}\right)^{-1}, \] and \[ m_n=s_n^2\left(\frac{m_0}{s_0^2}+\frac{n\bar y_n}{\sigma^2}\right). \]
Thus, posterior mean is a precision-weighted average of prior mean and sample mean.
set.seed(2026)
mu_true <- 2
sigma <- 1
nmax <- 600
y <- rnorm(nmax, mu_true, sigma)
n_grid2 <- c(1, 2, 5, 10, 20, 50, 100, 300, 600)
priors <- data.frame(prior = c("Vague prior", "Wrong but weak prior", "Wrong dogmatic prior"),
m0 = c(0, -3, -3),
s0 = c(10, 3, 0.25))
norm_tab <- data.frame()
for(i in seq_len(nrow(priors))) {
for(n in n_grid2) {
ybar <- mean(y[1:n])
s0 <- priors$s0[i]
m0 <- priors$m0[i]
sn2 <- 1/(1/s0^2 + n/sigma^2)
mn <- sn2 * (m0/s0^2 + n*ybar/sigma^2)
ci <- mn + c(-1,1)*qnorm(0.975)*sqrt(sn2)
norm_tab <- rbind(norm_tab, data.frame(prior = priors$prior[i], n=n,
post_mean=mn, ci_low=ci[1], ci_high=ci[2],
post_sd=sqrt(sn2)))
}
}
pretty_table(head(norm_tab, 12), digits = 4,
caption = "Posterior summaries for different priors as n increases.")| prior | n | post_mean | ci_low | ci_high | post_sd |
|---|---|---|---|---|---|
| Vague prior | 1 | 2.4956 | 0.5454 | 4.4459 | 0.9950 |
| Vague prior | 2 | 1.7119 | 0.3294 | 3.0943 | 0.7053 |
| Vague prior | 5 | 1.7622 | 0.8866 | 2.6379 | 0.4468 |
| Vague prior | 10 | 1.4183 | 0.7988 | 2.0378 | 0.3161 |
| Vague prior | 20 | 1.6057 | 1.1676 | 2.0439 | 0.2236 |
| Vague prior | 50 | 1.9833 | 1.7062 | 2.2605 | 0.1414 |
| Vague prior | 100 | 1.9018 | 1.7058 | 2.0978 | 0.1000 |
| Vague prior | 300 | 2.0416 | 1.9284 | 2.1547 | 0.0577 |
| Vague prior | 600 | 2.0567 | 1.9767 | 2.1367 | 0.0408 |
| Wrong but weak prior | 1 | 1.9685 | 0.1091 | 3.8279 | 0.9487 |
| Wrong but weak prior | 2 | 1.4720 | 0.1231 | 2.8209 | 0.6882 |
| Wrong but weak prior | 5 | 1.6621 | 0.7952 | 2.5291 | 0.4423 |
ggplot(norm_tab, aes(n, post_mean, color = prior, fill = prior)) +
geom_ribbon(aes(ymin = ci_low, ymax = ci_high), alpha = 0.10, color = NA) +
geom_line(linewidth = 1.1) + geom_point(size = 2) +
geom_hline(yintercept = mu_true, linetype = 2) +
scale_x_log10(breaks = n_grid2) +
labs(title = "Prior influence fades as data accumulate, unless the prior is very dogmatic",
subtitle = "Dashed line is the true mean used in simulation.",
x = "Sample size n, log scale", y = expression(posterior~mean~of~mu), color = "Prior", fill = "Prior")
6 6. Bayesian hypothesis learning: evidence accumulates
Bayesian hypothesis testing compares hypotheses through posterior probabilities or Bayes factors.
For a detector probability:
\[ H_0:\theta=0.5, \] versus \[ H_1:\theta\sim\text{Beta}(1,1). \]
The Bayes factor in favor of \(H_1\) is \[ BF_{10}=\frac{m_1(y)}{m_0(y)}, \] where \(m_1(y)\) and \(m_0(y)\) are marginal likelihoods under the two hypotheses.
set.seed(121)
theta_alt <- 0.58
nmax_bf <- 600
z <- rbinom(nmax_bf, 1, theta_alt)
n_bf_grid <- seq(10, nmax_bf, by = 10)
bf_tab <- data.frame()
for(n in n_bf_grid) {
yy <- sum(z[1:n])
# H1: Beta(1,1) prior, marginal likelihood is beta-binomial part ignoring combinatorial term.
logm1 <- lbeta(1 + yy, 1 + n - yy) - lbeta(1, 1)
# H0: theta = 0.5, likelihood part ignoring same combinatorial term.
logm0 <- yy * log(0.5) + (n - yy) * log(0.5)
logBF10 <- logm1 - logm0
bf_tab <- rbind(bf_tab, data.frame(n=n, y=yy, sample_prop=yy/n,
log10_BF10=logBF10/log(10),
post_prob_gt_half = 1 - pbeta(0.5, 1+yy, 1+n-yy)))
}
pretty_table(tail(bf_tab, 8), digits = 4,
caption = "Sequential Bayes factor against theta = 0.5 when true theta is 0.58.")| n | y | sample_prop | log10_BF10 | post_prob_gt_half | |
|---|---|---|---|---|---|
| 53 | 530 | 303 | 0.5717 | 1.1055 | 0.9995 |
| 54 | 540 | 307 | 0.5685 | 0.9361 | 0.9993 |
| 55 | 550 | 313 | 0.5691 | 1.0109 | 0.9994 |
| 56 | 560 | 318 | 0.5679 | 0.9660 | 0.9993 |
| 57 | 570 | 322 | 0.5649 | 0.8079 | 0.9990 |
| 58 | 580 | 328 | 0.5655 | 0.8807 | 0.9992 |
| 59 | 590 | 336 | 0.5695 | 1.1906 | 0.9996 |
| 60 | 600 | 343 | 0.5717 | 1.3899 | 0.9998 |
ggplot(bf_tab, aes(n, log10_BF10)) +
geom_line(linewidth = 1.1, color = "navy") +
geom_hline(yintercept = 0, linetype = 2) +
labs(title = "Bayesian evidence accumulates sequentially",
subtitle = "Positive log10(BF10) favors the flexible alternative over theta = 0.5.",
x = "Sample size n", y = expression(log[10](BF[10])))
ggplot(bf_tab, aes(n, post_prob_gt_half)) +
geom_line(linewidth = 1.1, color = "darkgreen") +
geom_hline(yintercept = 0.95, linetype = 2) +
labs(title = "Posterior probability of a positive departure from 0.5",
x = "Sample size n", y = expression(P(theta>0.5~"| data")))
7 7. Warning: larger data reveal truth only when the model can express it
A central lesson for scientific data science:
If the model is correct or close enough, the posterior can concentrate near the true parameter. If the model is misspecified, the posterior may concentrate near a pseudo-truth: the best approximation inside the wrong model class.
7.1 7.1 Misspecification simulation
Generate data from a mixture: \[ Y\sim 0.8N(0,1)+0.2N(5,1). \]
But fit a single Normal mean model. Then increasing data do not recover the dominant component mean 0; they estimate the mixture mean around 1.
set.seed(333)
Nmis <- 5000
component <- rbinom(Nmis, 1, 0.2)
y_mis <- rnorm(Nmis, mean = ifelse(component == 1, 5, 0), sd = 1)
mixture_mean <- mean(y_mis)
n_mis <- c(30, 100, 500, 2000, 5000)
mis_tab <- data.frame()
xmu <- seq(-1, 2.5, length.out = 500)
mis_dens <- list()
m0 <- 0; s0 <- 5; sig <- 1
for(n in n_mis) {
ybar <- mean(y_mis[1:n])
sn2 <- 1/(1/s0^2 + n/sig^2)
mn <- sn2*(m0/s0^2 + n*ybar/sig^2)
mis_dens[[as.character(n)]] <- data.frame(n = factor(n, levels = n_mis),
mu = xmu,
density = dnorm(xmu, mn, sqrt(sn2)))
ci <- mn + c(-1,1)*qnorm(0.975)*sqrt(sn2)
mis_tab <- rbind(mis_tab, data.frame(n=n, post_mean=mn, ci_low=ci[1], ci_high=ci[2]))
}
mis_dens <- do.call(rbind, mis_dens)
pretty_table(mis_tab, digits = 4,
caption = "Posterior under a wrong single-Normal model for mixture-generated data.")| n | post_mean | ci_low | ci_high |
|---|---|---|---|
| 30 | 0.9027 | 0.5451 | 1.2603 |
| 100 | 1.0974 | 0.9015 | 1.2934 |
| 500 | 0.9977 | 0.9100 | 1.0853 |
| 2000 | 1.0557 | 1.0118 | 1.0995 |
| 5000 | 1.0040 | 0.9763 | 1.0317 |
ggplot(data.frame(y = y_mis[1:1000]), aes(y)) +
geom_histogram(aes(y = after_stat(density)), bins = 45, fill = "grey80", color = "white") +
geom_density(linewidth = 1.1, color = "firebrick") +
geom_vline(xintercept = 0, linetype = 2, color = "black") +
geom_vline(xintercept = mean(y_mis), linetype = 3, color = "blue") +
labs(title = "The true data distribution is not a single Normal population",
subtitle = "Dashed black line: dominant component mean. Dotted blue line: mixture mean.",
x = "Observation", y = "Density")
ggplot(mis_dens, aes(mu, density, color = n)) +
geom_line(linewidth = 1.1) +
geom_vline(xintercept = 0, linetype = 2, color = "black") +
geom_vline(xintercept = mean(y_mis), linetype = 3, color = "blue") +
labs(title = "Under misspecification, posterior certainty can increase around the wrong scientific summary",
subtitle = "Posterior concentrates near mixture mean, not the dominant component mean.",
x = expression(mu), y = "Posterior density", color = "n")
This is one of the most important messages for research students: more data improve inference only relative to the model class one is willing to entertain.
8 8. Real data example I: Physical chemistry and vapor pressure
The R dataset pressure contains 19 observations of
temperature in degrees Celsius and vapor pressure of mercury in
millimeters of mercury. A scientific model should respect that vapor
pressure changes nonlinearly with temperature.
A common physical insight is that log pressure often becomes approximately linear in inverse absolute temperature.
data(pressure)
dfp <- pressure
dfp$tempK <- dfp$temperature + 273.15
dfp$invT1000 <- 1000/dfp$tempK
dfp$log_pressure <- log(dfp$pressure)
pretty_table(dfp, digits = 4, caption = "Mercury vapor pressure data.")| temperature | pressure | tempK | invT1000 | log_pressure |
|---|---|---|---|---|
| 0 | 0.0002 | 273.15 | 3.6610 | -8.5172 |
| 20 | 0.0012 | 293.15 | 3.4112 | -6.7254 |
| 40 | 0.0060 | 313.15 | 3.1934 | -5.1160 |
| 60 | 0.0300 | 333.15 | 3.0017 | -3.5066 |
| 80 | 0.0900 | 353.15 | 2.8317 | -2.4079 |
| 100 | 0.2700 | 373.15 | 2.6799 | -1.3093 |
| 120 | 0.7500 | 393.15 | 2.5436 | -0.2877 |
| 140 | 1.8500 | 413.15 | 2.4204 | 0.6152 |
| 160 | 4.2000 | 433.15 | 2.3087 | 1.4351 |
| 180 | 8.8000 | 453.15 | 2.2068 | 2.1748 |
| 200 | 17.3000 | 473.15 | 2.1135 | 2.8507 |
| 220 | 32.1000 | 493.15 | 2.0278 | 3.4689 |
| 240 | 57.0000 | 513.15 | 1.9487 | 4.0431 |
| 260 | 96.0000 | 533.15 | 1.8756 | 4.5643 |
| 280 | 157.0000 | 553.15 | 1.8078 | 5.0562 |
| 300 | 247.0000 | 573.15 | 1.7447 | 5.5094 |
| 320 | 376.0000 | 593.15 | 1.6859 | 5.9296 |
| 340 | 558.0000 | 613.15 | 1.6309 | 6.3244 |
| 360 | 806.0000 | 633.15 | 1.5794 | 6.6921 |
p1 <- ggplot(dfp, aes(temperature, pressure)) +
geom_point(size = 3, color = "firebrick") +
labs(title = "Raw scale: strong nonlinearity",
x = "Temperature (deg C)", y = "Pressure (mm Hg)")
p2 <- ggplot(dfp, aes(invT1000, log_pressure)) +
geom_point(size = 3, color = "navy") +
labs(title = "Physics-aware transform: log pressure vs inverse temperature",
x = expression(1000/T[K]), y = "log pressure")
p1
8.1 8.1 Bayesian linear model on transformed scale
Fit \[ \log(P_i)=\beta_0+\beta_1\frac{1000}{T_i}+\epsilon_i, \qquad \epsilon_i\sim N(0,\sigma^2). \]
X <- cbind(1, dfp$invT1000)
y <- dfp$log_pressure
fit_p <- bayes_lm_nig(X, y, a0 = 2, b0 = 1)
samp_p <- sample_bayes_lm(fit_p, 4000)
colnames(samp_p$beta) <- c("intercept", "invT1000")
coef_tab <- data.frame(parameter = colnames(samp_p$beta),
t(apply(samp_p$beta, 2, qfun)))
pretty_table(coef_tab, digits = 4, caption = "Posterior coefficient summaries for transformed vapor-pressure model.")| parameter | mean | q025 | q50 | q975 | |
|---|---|---|---|---|---|
| intercept | intercept | 18.2809 | 17.7251 | 18.2781 | 18.8328 |
| invT1000 | invT1000 | -7.3090 | -7.5399 | -7.3093 | -7.0802 |
gridp <- data.frame(invT1000 = seq(min(dfp$invT1000), max(dfp$invT1000), length.out = 120))
Xg <- cbind(1, gridp$invT1000)
mu_draw <- Xg %*% t(samp_p$beta)
gridp$mean <- rowMeans(mu_draw)
gridp$lo <- apply(mu_draw, 1, quantile, 0.025)
gridp$hi <- apply(mu_draw, 1, quantile, 0.975)
ggplot(dfp, aes(invT1000, log_pressure)) +
geom_ribbon(data = gridp, aes(x = invT1000, ymin = lo, ymax = hi), alpha = 0.18, fill = "darkgreen", inherit.aes = FALSE) +
geom_line(data = gridp, aes(x = invT1000, y = mean), linewidth = 1.1, color = "darkgreen", inherit.aes = FALSE) +
geom_point(size = 3, color = "navy") +
labs(title = "Bayesian posterior mean and credible band for vapor pressure law",
subtitle = "A physically motivated transform gives a simple and interpretable statistical model.",
x = expression(1000/T[K]), y = "log pressure")
9 9. Real data example II: Nile flow changepoint as Bayesian model uncertainty
The Nile dataset contains annual flow at Aswan from 1871
to 1970. It is famous for an apparent changepoint near 1898.
We model an unknown changepoint \(\tau\). Before \(\tau\), the mean flow is one level; after \(\tau\), the mean flow is another level.
yN <- as.numeric(Nile)
yearsN <- as.numeric(time(Nile))
dfn <- data.frame(year = yearsN, flow = yN)
ggplot(dfn, aes(year, flow)) +
geom_line(linewidth = 0.8, color = "steelblue") +
geom_point(size = 1.7, color = "steelblue") +
labs(title = "Nile annual flow: hydrological time series with possible level shift",
x = "Year", y = expression(flow~(10^8~m^3)))
9.1 9.1 Posterior distribution over changepoint
We use a simple Bayesian grid calculation: uniform prior over possible changepoints and integrated Normal likelihood for the two segments.
seg_log_marginal <- function(v) {
n <- length(v)
if(n < 3) return(-Inf)
sse <- sum((v - mean(v))^2)
sse <- max(sse, 1e-8)
# Proportional integrated likelihood for Normal data with unknown mean and variance.
lgamma((n-1)/2) - 0.5*log(n) - ((n-1)/2)*log(pi*sse)
}
candidates <- yearsN[8:(length(yearsN)-8)]
logpost <- numeric(length(candidates))
for(i in seq_along(candidates)) {
tau <- candidates[i]
before <- yN[yearsN <= tau]
after <- yN[yearsN > tau]
logpost[i] <- seg_log_marginal(before) + seg_log_marginal(after)
}
post_tau <- exp(logpost - max(logpost))
post_tau <- post_tau / sum(post_tau)
tau_df <- data.frame(tau = candidates, post = post_tau)
tau_map <- tau_df$tau[which.max(tau_df$post)]
tau_mean <- sum(tau_df$tau * tau_df$post)
pretty_table(data.frame(MAP_changepoint = tau_map,
posterior_mean_changepoint = tau_mean),
digits = 2, caption = "Posterior summary for Nile changepoint.")| MAP_changepoint | posterior_mean_changepoint |
|---|---|
| 1898 | 1897.86 |
ggplot(tau_df, aes(tau, post)) +
geom_col(fill = "darkorange", alpha = 0.8) +
geom_vline(xintercept = tau_map, linetype = 2) +
labs(title = "Posterior distribution over possible Nile changepoint",
subtitle = "The posterior assigns probability to many plausible changepoints, not only one estimate.",
x = "Candidate changepoint year", y = "Posterior probability")
dfn$period <- ifelse(dfn$year <= tau_map, "Before MAP changepoint", "After MAP changepoint")
means_nile <- aggregate(flow ~ period, dfn, mean)
ggplot(dfn, aes(year, flow)) +
geom_line(color = "grey55") +
geom_point(aes(color = period), size = 2) +
geom_vline(xintercept = tau_map, linetype = 2) +
geom_hline(data = means_nile, aes(yintercept = flow, color = period), linewidth = 1.1) +
labs(title = "Nile flow with posterior MAP changepoint",
subtitle = paste("MAP changepoint:", tau_map),
x = "Year", y = expression(flow~(10^8~m^3)), color = "")
This example is conceptually important: the Bayesian output is not merely a detected year. It is a posterior distribution over possible changepoint years.
10 10. Real data example III: Solar activity and Bayesian cyclic regression
The sunspot.year dataset contains yearly sunspot numbers
from 1700 to 1988.
We fit an intentionally simple cyclic regression to log-transformed sunspot counts: \[ \log(1+Y_t)=\beta_0+\beta_1\sin(2\pi t/11)+\beta_2\cos(2\pi t/11)+\epsilon_t. \]
This is not a final solar-physics model. It is an inference demonstration: scientific structure can enter through features.
dfs <- data.frame(year = as.numeric(time(sunspot.year)),
sunspot = as.numeric(sunspot.year))
dfs$t <- dfs$year - min(dfs$year)
dfs$log_sunspot <- log1p(dfs$sunspot)
dfs$sin11 <- sin(2*pi*dfs$t/11)
dfs$cos11 <- cos(2*pi*dfs$t/11)
dfs$sin22 <- sin(2*pi*dfs$t/22)
dfs$cos22 <- cos(2*pi*dfs$t/22)
pretty_table(head(dfs, 10), digits = 3, caption = "First rows of yearly sunspot data with cyclic features.")| year | sunspot | t | log_sunspot | sin11 | cos11 | sin22 | cos22 |
|---|---|---|---|---|---|---|---|
| 1700 | 5 | 0 | 1.792 | 0.000 | 1.000 | 0.000 | 1.000 |
| 1701 | 11 | 1 | 2.485 | 0.541 | 0.841 | 0.282 | 0.959 |
| 1702 | 16 | 2 | 2.833 | 0.910 | 0.415 | 0.541 | 0.841 |
| 1703 | 23 | 3 | 3.178 | 0.990 | -0.142 | 0.756 | 0.655 |
| 1704 | 36 | 4 | 3.611 | 0.756 | -0.655 | 0.910 | 0.415 |
| 1705 | 58 | 5 | 4.078 | 0.282 | -0.959 | 0.990 | 0.142 |
| 1706 | 29 | 6 | 3.401 | -0.282 | -0.959 | 0.990 | -0.142 |
| 1707 | 20 | 7 | 3.045 | -0.756 | -0.655 | 0.910 | -0.415 |
| 1708 | 10 | 8 | 2.398 | -0.990 | -0.142 | 0.756 | -0.655 |
| 1709 | 8 | 9 | 2.197 | -0.910 | 0.415 | 0.541 | -0.841 |
ggplot(dfs, aes(year, sunspot)) +
geom_line(color = "firebrick", linewidth = 0.7) +
labs(title = "Yearly sunspot numbers",
x = "Year", y = "Sunspot number")
Xs <- as.matrix(cbind(1, dfs$sin11, dfs$cos11, dfs$sin22, dfs$cos22))
ys <- dfs$log_sunspot
fit_s <- bayes_lm_nig(Xs, ys, a0 = 2, b0 = 1)
samp_s <- sample_bayes_lm(fit_s, 3000)
colnames(samp_s$beta) <- c("intercept", "sin11", "cos11", "sin22", "cos22")
sun_coef <- data.frame(parameter = colnames(samp_s$beta),
t(apply(samp_s$beta, 2, qfun)))
pretty_table(sun_coef, digits = 4, caption = "Posterior coefficient summaries for cyclic sunspot regression.")| parameter | mean | q025 | q50 | q975 | |
|---|---|---|---|---|---|
| intercept | intercept | 3.5066 | 3.4090 | 3.5070 | 3.6038 |
| sin11 | sin11 | -0.3012 | -0.4304 | -0.2991 | -0.1707 |
| cos11 | cos11 | -0.7633 | -0.8985 | -0.7636 | -0.6256 |
| sin22 | sin22 | -0.0474 | -0.1877 | -0.0460 | 0.0946 |
| cos22 | cos22 | 0.1145 | -0.0224 | 0.1135 | 0.2523 |
mu_s <- Xs %*% t(samp_s$beta)
dfs$fit <- rowMeans(mu_s)
dfs$lo <- apply(mu_s, 1, quantile, 0.025)
dfs$hi <- apply(mu_s, 1, quantile, 0.975)
ggplot(dfs, aes(year, log_sunspot)) +
geom_ribbon(aes(ymin = lo, ymax = hi), alpha = 0.18, fill = "steelblue") +
geom_line(aes(y = fit), color = "navy", linewidth = 1) +
geom_line(color = "grey40", alpha = 0.7) +
labs(title = "Bayesian cyclic regression for solar activity",
subtitle = "The model is simple but demonstrates posterior uncertainty around a scientific feature representation.",
x = "Year", y = "log(1 + sunspot number)")
11 11. Real data example IV: Earthquake rare-event probability with Bayesian uncertainty wrapper
The quakes dataset contains 1000 seismic events near
Fiji, with latitude, longitude, depth, magnitude, and number of
reporting stations.
We define a toy scientific task: estimate the probability that an event is relatively large, here \(\text{magnitude}\ge 5.0\), using depth and reporting stations.
This is not a hazard model for policy use; it is an inference demonstration.
data(quakes)
dfq <- quakes
dfq$large <- as.integer(dfq$mag >= 5.0)
dfq$depth_z <- as.numeric(scale(dfq$depth))
dfq$stations_z <- as.numeric(scale(dfq$stations))
pretty_table(head(dfq, 10), digits = 3, caption = "First rows of Fiji earthquake data.")| lat | long | depth | mag | stations | large | depth_z | stations_z |
|---|---|---|---|---|---|---|---|
| -20.42 | 181.62 | 562 | 4.8 | 41 | 0 | 1.163 | 0.346 |
| -20.62 | 181.03 | 650 | 4.2 | 15 | 0 | 1.571 | -0.841 |
| -26.00 | 184.10 | 42 | 5.4 | 43 | 1 | -1.250 | 0.438 |
| -17.97 | 181.66 | 626 | 4.1 | 19 | 0 | 1.460 | -0.658 |
| -20.42 | 181.96 | 649 | 4.0 | 11 | 0 | 1.566 | -1.024 |
| -19.68 | 184.31 | 195 | 4.0 | 12 | 0 | -0.540 | -0.978 |
| -11.70 | 166.10 | 82 | 4.8 | 43 | 0 | -1.064 | 0.438 |
| -28.11 | 181.93 | 194 | 4.4 | 15 | 0 | -0.545 | -0.841 |
| -28.74 | 181.74 | 211 | 4.7 | 35 | 0 | -0.466 | 0.072 |
| -17.47 | 179.59 | 622 | 4.3 | 19 | 0 | 1.441 | -0.658 |
ggplot(dfq, aes(depth, mag)) +
geom_point(aes(color = factor(large)), alpha = 0.65, size = 1.9) +
labs(title = "Fiji earthquakes: magnitude versus depth",
subtitle = "Large event indicator defined here as magnitude >= 5.0 for demonstration.",
x = "Depth (km)", y = "Magnitude", color = "Large")
fit_glm <- glm(large ~ depth_z + stations_z, data = dfq, family = binomial())
bhat <- coef(fit_glm)
Vhat <- vcov(fit_glm)
set.seed(77)
bdraw <- MASS::mvrnorm(4000, mu = bhat, Sigma = Vhat)
colnames(bdraw) <- names(bhat)
logit_tab <- data.frame(parameter = colnames(bdraw), t(apply(bdraw, 2, qfun)))
pretty_table(logit_tab, digits = 4,
caption = "Approximate Bayesian posterior summaries using large-sample Normal approximation to logistic regression.")| parameter | mean | q025 | q50 | q975 | |
|---|---|---|---|---|---|
| (Intercept) | (Intercept) | -2.6189 | -2.9626 | -2.6195 | -2.2643 |
| depth_z | depth_z | -0.4001 | -0.6587 | -0.4032 | -0.1274 |
| stations_z | stations_z | 3.1748 | 2.7297 | 3.1734 | 3.6216 |
gridq <- data.frame(depth = seq(min(dfq$depth), max(dfq$depth), length.out = 160))
gridq$depth_z <- (gridq$depth - mean(dfq$depth))/sd(dfq$depth)
gridq$stations_z <- 0
Xq <- model.matrix(~ depth_z + stations_z, data = gridq)
pdraw <- inv_logit(Xq %*% t(bdraw))
gridq$pmean <- rowMeans(pdraw)
gridq$plo <- apply(pdraw, 1, quantile, 0.025)
gridq$phi <- apply(pdraw, 1, quantile, 0.975)
ggplot(gridq, aes(depth, pmean)) +
geom_ribbon(aes(ymin = plo, ymax = phi), alpha = 0.18, fill = "darkgreen") +
geom_line(linewidth = 1.1, color = "darkgreen") +
labs(title = "Posterior uncertainty for probability of a relatively large earthquake",
subtitle = "Stations fixed at average level; uncertainty band comes from approximate posterior draws.",
x = "Depth (km)", y = expression(P(magnitude >= 5.0)))
This example links ML-style classification and statistical inference. A classifier can produce probabilities. Bayesian or approximate Bayesian inference adds uncertainty around those probabilities.
12 12. Real data example V: Galaxy velocities and model uncertainty
The MASS::galaxies dataset contains velocities of 82
galaxies in the Corona Borealis region. Multimodality in such surveys is
scientifically relevant because it can indicate voids and
superclusters.
We fit Normal mixtures with \(K=1,2,3,4\) components using EM, then use BIC weights as a rough model-uncertainty summary.
normal_mix_em <- function(x, K, nstart = 8, maxit = 300, tol = 1e-7) {
x <- as.numeric(x)
n <- length(x)
best <- NULL
for(st in seq_len(nstart)) {
if(K == 1) {
pi <- 1
mu <- mean(x)
sig <- sd(x)
} else {
mu <- as.numeric(quantile(x, probs = seq(0.15, 0.85, length.out = K))) + rnorm(K, 0, sd(x)/20)
sig <- rep(sd(x)/K, K)
pi <- rep(1/K, K)
}
oldll <- -Inf
for(iter in seq_len(maxit)) {
dens <- sapply(seq_len(K), function(k) pi[k] * dnorm(x, mu[k], sig[k]))
if(K == 1) dens <- matrix(dens, ncol = 1)
denom <- rowSums(dens) + 1e-300
resp <- dens / denom
Nk <- colSums(resp)
pi <- as.numeric(Nk/n)
mu <- as.numeric(colSums(resp * x)/Nk)
sig <- sqrt(as.numeric(colSums(resp * (x - rep(mu, each = n))^2)/Nk))
sig <- pmax(sig, sd(x)/100)
ll <- sum(log(denom))
if(abs(ll - oldll) < tol) break
oldll <- ll
}
fit <- list(K=K, pi=pi, mu=mu, sig=sig, logLik=ll, iter=iter)
if(is.null(best) || fit$logLik > best$logLik) best <- fit
}
best
}gal <- as.numeric(MASS::galaxies)/1000
fit_list <- lapply(1:4, function(K) normal_mix_em(gal, K))
bic_tab <- data.frame()
for(K in 1:4) {
fit <- fit_list[[K]]
df <- (K-1) + K + K
bic <- -2*fit$logLik + df*log(length(gal))
bic_tab <- rbind(bic_tab, data.frame(K=K, logLik=fit$logLik, df=df, BIC=bic))
}
bic_tab$approx_model_weight <- exp(-0.5*(bic_tab$BIC - min(bic_tab$BIC)))
bic_tab$approx_model_weight <- bic_tab$approx_model_weight/sum(bic_tab$approx_model_weight)
pretty_table(bic_tab, digits = 4,
caption = "Mixture model comparison for galaxy velocities using BIC weights.")| K | logLik | df | BIC | approx_model_weight |
|---|---|---|---|---|
| 1 | -240.3379 | 2 | 489.4892 | 0.0000 |
| 2 | -220.2433 | 5 | 462.5202 | 0.0000 |
| 3 | -203.1792 | 8 | 441.6122 | 0.9360 |
| 4 | -199.2527 | 11 | 446.9793 | 0.0639 |
xg <- seq(min(gal)-2, max(gal)+2, length.out = 600)
dens_mix <- data.frame()
for(K in 1:4) {
fit <- fit_list[[K]]
dd <- rep(0, length(xg))
for(k in seq_len(K)) dd <- dd + fit$pi[k]*dnorm(xg, fit$mu[k], fit$sig[k])
dens_mix <- rbind(dens_mix, data.frame(K = factor(K), velocity = xg, density = dd))
}
ggplot(data.frame(velocity = gal), aes(velocity)) +
geom_histogram(aes(y = after_stat(density)), bins = 24, fill = "grey85", color = "white") +
geom_line(data = dens_mix, aes(velocity, density, color = K), linewidth = 1.1) +
labs(title = "Galaxy velocities: mixture models express multimodal scientific structure",
subtitle = "BIC weights approximate model uncertainty over number of velocity groups.",
x = "Velocity (1000 km/s)", y = "Density", color = "K components")
This demonstrates an important Bayesian-data-science principle: sometimes the scientific question is not only a parameter value but also model structure.
13 13. Posterior predictive distribution: predicting future data
Bayesian inference naturally predicts future observations by averaging over posterior uncertainty: \[ p(\widetilde y\mid y)=\int p(\widetilde y\mid\theta)\pi(\theta\mid y)d\theta. \]
In the vapor pressure example, posterior predictive uncertainty combines uncertainty in regression coefficients and residual noise.
gridp2 <- data.frame(invT1000 = seq(min(dfp$invT1000), max(dfp$invT1000), length.out = 100))
Xg2 <- cbind(1, gridp2$invT1000)
ndraw <- 2500
idx <- sample(seq_len(nrow(samp_p$beta)), ndraw)
mu_pred <- Xg2 %*% t(samp_p$beta[idx, ])
# Include observation-level noise for posterior predictive draws
sigma_draw <- sqrt(samp_p$sigma2[idx])
yrep <- mu_pred
for(j in seq_len(ncol(yrep))) {
yrep[, j] <- yrep[, j] + rnorm(nrow(yrep), 0, sigma_draw[j])
}
gridp2$pred_mean <- rowMeans(yrep)
gridp2$pred_lo <- apply(yrep, 1, quantile, 0.025)
gridp2$pred_hi <- apply(yrep, 1, quantile, 0.975)
ggplot(dfp, aes(invT1000, log_pressure)) +
geom_ribbon(data = gridp2, aes(x = invT1000, ymin = pred_lo, ymax = pred_hi),
inherit.aes = FALSE, alpha = 0.18, fill = "purple") +
geom_line(data = gridp2, aes(x = invT1000, y = pred_mean), color = "purple", linewidth = 1.1, inherit.aes = FALSE) +
geom_point(size = 3, color = "black") +
labs(title = "Posterior predictive band: uncertainty for future measurements",
subtitle = "Wider than the mean credible band because it includes measurement noise.",
x = expression(1000/T[K]), y = "log pressure")
14 14. Bayesian workflow for basic science
A mature Bayesian analysis is not only prior times likelihood. It is a workflow:
- Define the scientific estimand.
- Specify a probability model for data and unknowns.
- Choose priors that are scientifically meaningful or weakly informative.
- Compute the posterior.
- Verify with simulation when possible.
- Check posterior predictive adequacy.
- Test sensitivity to priors and likelihood choices.
- Report uncertainty clearly.
- Translate posterior summaries into scientific conclusions.
15 15. Summary
This RPubs note showed how general statistical inference connects to Bayesian inference.
Key lessons:
- Bayesian inference treats unknown scientific quantities probabilistically.
- The posterior is a conditional probability distribution over parameters after observing data.
- In correctly specified models, posterior uncertainty shrinks and posterior mass tends to concentrate near the truth as data increase.
- Priors matter most when data are scarce; regular priors usually become less influential as information grows.
- Under misspecification, larger data can make the wrong model confidently wrong.
- Real science needs uncertainty-aware models, model checking, and scientific structure.
16 16. Suggested next topics
The next natural topics are:
- Bayesian computation: MCMC, Hamiltonian Monte Carlo, variational Bayes.
- Bayesian hierarchical models.
- Bayesian nonparametrics.
- Bayesian model comparison and Bayes factors.
- Bayesian decision theory.
- Bayesian machine learning.
- Bayesian inverse problems in astronomy, geoscience, soft matter, and biological sciences.
17 References and data notes
The examples use standard R datasets and MASS::galaxies.
The pressure dataset contains 19 measurements of mercury
vapor pressure against temperature. The Nile dataset
contains annual Nile flow at Aswan from 1871–1970 and is known for an
apparent changepoint near 1898. The sunspot.year dataset
contains yearly sunspot numbers from 1700–1988. The quakes
dataset contains 1000 seismic events near Fiji with magnitude greater
than 4.0. The MASS::galaxies dataset contains velocities of
82 galaxies in the Corona Borealis region.