1 Aim of this Part II

This note is a continuation of a first lecture note on Bayesian statistical inference. The goal here is to move from the basic slogan

\[ \text{posterior} \propto \text{likelihood}\times \text{prior} \]

to the deeper research-level workflow used in scientific data analysis:

  1. define the probability environment carefully;
  2. specify the sampling model and prior;
  3. compute the posterior analytically when possible;
  4. compute the posterior approximately when necessary;
  5. check the model by posterior prediction;
  6. compare models without pretending that any model is literally true;
  7. use ML for structure discovery and Bayesian inference for uncertainty.

The examples use astronomy-flavoured and basic-science situations: photon counting, sunspot time series, galaxy velocities, robust Hubble-type regression and hierarchical shrinkage of repeated scientific measurements.

2 The probability environment behind Bayesian inference

2.1 Data space, parameter space, model and prior

Let the observable data be a random object

\[ Y:(\Omega,\mathcal F,P)\longrightarrow (\mathcal Y,\mathcal A), \]

where:

  • \(\Omega\) is the underlying sample space of possible worlds;
  • \(\mathcal F\) is the sigma-field of events;
  • \(P\) is the probability law;
  • \((\mathcal Y,\mathcal A)\) is the measurable data space.

A parametric statistical model is a family of probability measures

\[ \mathcal P=\{P_\theta:\theta\in\Theta\} \]

on \((\mathcal Y,\mathcal A)\). If the family is dominated by a reference measure \(\nu\), then

\[ p(y\mid \theta)=\frac{dP_\theta}{d\nu}(y) \]

is the likelihood density or mass function as a function of \(y\) and \(\theta\).

Bayesian inference adds a prior probability measure \(\Pi\) on the parameter space \((\Theta,\mathcal T)\). The joint law is

\[ P(dy,d\theta)=P_\theta(dy)\Pi(d\theta). \]

The posterior is the conditional distribution of \(\theta\) given the observed data \(Y=y\):

\[ \Pi(B\mid y)= \frac{\int_B p(y\mid \theta)\,\Pi(d\theta)} {\int_\Theta p(y\mid \theta)\,\Pi(d\theta)}, \qquad B\in\mathcal T. \]

The denominator

\[ m(y)=\int_\Theta p(y\mid \theta)\,\Pi(d\theta) \]

is the marginal likelihood or evidence. It is essential for Bayes factors and model comparison.

2.2 Posterior predictive distribution

Scientific inference often asks not only “what is the parameter?” but also “what future or unseen data should we expect?” The posterior predictive distribution is

\[ p(y_{\rm new}\mid y)=\int_\Theta p(y_{\rm new}\mid \theta)\,\Pi(d\theta\mid y). \]

This is the Bayesian mechanism for propagating parameter uncertainty into prediction.

2.3 Bayesian decision rule

If an action \(a\in\mathcal A_0\) is chosen under loss \(L(a,\theta)\), the Bayes action is

\[ a^\star(y)=\arg\min_{a\in\mathcal A_0} E\{L(a,\theta)\mid y\}. \]

For squared-error loss, the posterior mean is the Bayes action. For absolute-error loss, the posterior median is the Bayes action. For 0–1 classification loss, the posterior mode/class with largest posterior probability is the Bayes action.

3 Exponential families and conjugacy

3.1 Why conjugacy matters

A conjugate prior gives a posterior in the same family as the prior. Conjugacy is not required for Bayesian inference, but it is extremely useful pedagogically and computationally. It shows how data update prior information in closed form.

A regular exponential family has density/mass

\[ p(y\mid \theta)=h(y)\exp\{\eta(\theta)^\top T(y)-A(\theta)\}, \]

where:

  • \(T(y)\) is the sufficient statistic;
  • \(\eta(\theta)\) is the natural parameter;
  • \(A(\theta)\) is the log-normalizing function.

Many conjugate priors work because the likelihood depends on data through low-dimensional sufficient statistics.

conj <- data.frame(
  Scientific_situation = c("Binary detection", "Photon / event counts", "Gaussian measurement with known variance", "Gaussian measurement with unknown variance", "Category probabilities"),
  Sampling_model = c("Binomial(n, theta)", "Poisson(lambda)", "Normal(mu, sigma^2 known)", "Normal(mu, sigma^2 unknown)", "Multinomial"),
  Conjugate_prior = c("Beta(a,b)", "Gamma(a,b)", "Normal(m0, v0)", "Normal-Inverse-Gamma", "Dirichlet(alpha)"),
  Posterior_update = c("a + successes, b + failures", "a + total count, b + exposure", "precision adds", "mean and variance update", "alpha_j + count_j"),
  Natural_science_use = c("instrument detection probability", "astronomical photon counts", "calibration measurements", "laboratory measurement uncertainty", "composition / class fractions")
)
knitr::kable(conj, caption = "Common conjugate Bayesian models in scientific inference.")
Common conjugate Bayesian models in scientific inference.
Scientific_situation Sampling_model Conjugate_prior Posterior_update Natural_science_use
Binary detection Binomial(n, theta) Beta(a,b) a + successes, b + failures instrument detection probability
Photon / event counts Poisson(lambda) Gamma(a,b) a + total count, b + exposure astronomical photon counts
Gaussian measurement with known variance Normal(mu, sigma^2 known) Normal(m0, v0) precision adds calibration measurements
Gaussian measurement with unknown variance Normal(mu, sigma^2 unknown) Normal-Inverse-Gamma mean and variance update laboratory measurement uncertainty
Category probabilities Multinomial Dirichlet(alpha) alpha_j + count_j composition / class fractions

4 Simulation 1: Bayesian updating is learning by reweighting uncertainty

Consider a detection experiment. A telescope pipeline detects a weak object with unknown probability \(\theta\). We observe \(y\) detections in \(n\) repeated trials. The model is

\[ Y\mid \theta\sim {\rm Binomial}(n,\theta), \qquad \theta\sim {\rm Beta}(a,b). \]

The posterior is

\[ \theta\mid Y=y\sim {\rm Beta}(a+y,b+n-y). \]

true_theta <- 0.27
n <- 60
y <- rbinom(1, n, true_theta)
a0 <- 2
b0 <- 8
apost <- a0 + y
bpost <- b0 + n - y

th <- seq(0, 1, length.out = 1000)
df_beta <- data.frame(
  theta = rep(th, 3),
  density = c(dbeta(th, a0, b0),
              dbeta(th, apost, bpost),
              dbeta(th, y + 1, n - y + 1)),
  curve = rep(c("Prior Beta(2,8)", "Posterior", "Likelihood shape"), each = length(th))
)

post_samp <- rbeta(20000, apost, bpost)
ci_beta <- quantile(post_samp, c(0.025, 0.5, 0.975))
summary_beta <- data.frame(
  n = n,
  detections = y,
  sample_fraction = y/n,
  posterior_mean = apost/(apost + bpost),
  posterior_median = ci_beta[2],
  lower_95 = ci_beta[1],
  upper_95 = ci_beta[3],
  true_theta = true_theta
)
knitr::kable(summary_beta, digits = 3, caption = "Beta--Binomial posterior summary.")
Beta–Binomial posterior summary.
n detections sample_fraction posterior_mean posterior_median lower_95 upper_95 true_theta
50% 60 13 0.217 0.214 0.212 0.129 0.317 0.27 
ggplot(df_beta, aes(theta, density, colour = curve)) +
  geom_line(linewidth = 1.1) +
  geom_vline(xintercept = true_theta, linetype = 2) +
  labs(title = "Bayesian learning in a Binomial detection experiment",
       subtitle = "Dashed vertical line = true detection probability used in simulation",
       x = expression(theta), y = "density", colour = "Curve")

Interpretation. The prior is not a fixed conclusion; it is an uncertainty distribution. The likelihood reshapes it using data. The posterior is a compromise, with the amount of compromise controlled by the amount of information in the data and prior.

5 Simulation 2: Photon counting with background contamination

5.1 Scientific motivation

In astronomy, we often observe photons in a source aperture. The count in the source aperture contains both source photons and background photons. A nearby background aperture estimates the background level.

Let:

\[ Y\mid s,b\sim {\rm Poisson}(t_s(s+b)), \qquad B\mid b\sim {\rm Poisson}(t_b b), \]

where:

  • \(s\geq 0\) is the source rate;
  • \(b\geq 0\) is the background rate;
  • \(t_s,t_b\) are exposure/area factors;
  • \(Y\) is the source-aperture count;
  • \(B\) is the background-aperture count.

A simple Bayesian model is

\[ s\sim {\rm Gamma}(a_s,r_s),\qquad b\sim {\rm Gamma}(a_b,r_b), \]

independently. Because \(s+b\) appears in the likelihood, the posterior of \((s,b)\) is not a simple product of gammas, but it can be computed accurately on a grid for teaching.

# True rates and exposures
s_true <- 4.0
b_true <- 1.7
ts <- 12
tb <- 40
Y <- rpois(1, ts * (s_true + b_true))
B <- rpois(1, tb * b_true)

# Independent gamma priors, shape-rate parameterization
as <- 1.5; rs <- 0.5
ab <- 1.5; rb <- 0.5

s_grid <- seq(0.001, 10, length.out = 260)
b_grid <- seq(0.001, 6, length.out = 240)

logpost <- outer(s_grid, b_grid, Vectorize(function(s, b) {
  dpois(Y, ts * (s + b), log = TRUE) +
    dpois(B, tb * b, log = TRUE) +
    dgamma(s, shape = as, rate = rs, log = TRUE) +
    dgamma(b, shape = ab, rate = rb, log = TRUE)
}))
post <- exp(logpost - max(logpost))
post <- post / sum(post)

marg_s <- rowSums(post); marg_s <- marg_s / sum(marg_s)
marg_b <- colSums(post); marg_b <- marg_b / sum(marg_b)

q_from_grid <- function(grid, prob, p) grid[which(cumsum(prob) >= p)[1]]
source_ci <- c(q_from_grid(s_grid, marg_s, 0.025),
               q_from_grid(s_grid, marg_s, 0.50),
               q_from_grid(s_grid, marg_s, 0.975))
back_ci <- c(q_from_grid(b_grid, marg_b, 0.025),
             q_from_grid(b_grid, marg_b, 0.50),
             q_from_grid(b_grid, marg_b, 0.975))

photon_summary <- data.frame(
  quantity = c("source count Y", "background count B", "source rate s", "background rate b"),
  observed_or_true = c(Y, B, s_true, b_true),
  posterior_median = c(NA, NA, source_ci[2], back_ci[2]),
  lower_95 = c(NA, NA, source_ci[1], back_ci[1]),
  upper_95 = c(NA, NA, source_ci[3], back_ci[3])
)
knitr::kable(photon_summary, digits = 3, caption = "Photon-counting posterior summary.")
Photon-counting posterior summary.
quantity observed_or_true posterior_median lower_95 upper_95
source count Y 70.0 NA NA NA
background count B 65.0 NA NA NA
source rate s 4.0 4.055 2.781 5.522
background rate b 1.7 1.658 1.281 2.084 
post_df <- expand.grid(source_rate = s_grid, background_rate = b_grid)
post_df$posterior_mass <- as.vector(post)

ggplot(post_df, aes(source_rate, background_rate, fill = posterior_mass)) +
  geom_raster(interpolate = TRUE) +
  geom_point(aes(x = s_true, y = b_true), colour = "white", size = 3) +
  scale_fill_viridis_c(option = "C") +
  labs(title = "Joint posterior over source and background rates",
       subtitle = "White point = true simulated value",
       x = "source rate s", y = "background rate b", fill = "posterior mass")

marg_df <- rbind(
  data.frame(value = s_grid, probability = marg_s, parameter = "source rate s"),
  data.frame(value = b_grid, probability = marg_b, parameter = "background rate b")
)

ggplot(marg_df, aes(value, probability)) +
  geom_line(linewidth = 1.1) +
  facet_wrap(~ parameter, scales = "free") +
  labs(title = "Marginal posterior distributions from the photon-counting model",
       x = "rate", y = "posterior probability on grid")

Research lesson. Background is not just a nuisance to be subtracted. It is an uncertain latent quantity. Bayesian inference lets us propagate uncertainty in background into uncertainty in source intensity.

6 Posterior predictive checking: a posterior can be internally valid but scientifically inadequate

6.1 Example: sunspot counts

The yearly sunspot.year data in R contains 289 yearly sunspot-number observations from 1700 to 1988. Monthly sunspot data are also available in R through sunspot.month, a longer and updated time series. We use the yearly series here for a simple reproducible example.

A naive model might say:

\[ Y_t\mid \lambda\sim {\rm Poisson}(\lambda),\qquad \lambda\sim {\rm Gamma}(a,b). \]

This model has a perfectly valid posterior, but it ignores solar cycles, time dependence and overdispersion.

y_sun <- as.numeric(sunspot.year)
y_sun_int <- round(y_sun)
n_sun <- length(y_sun_int)
year_sun <- as.numeric(time(sunspot.year))

# Gamma prior for Poisson rate
# Weakly informative around broad solar-count scale
a0 <- 2
b0 <- 0.02
apost <- a0 + sum(y_sun_int)
bpost <- b0 + n_sun

S <- 1500
lambda_draw <- rgamma(S, shape = apost, rate = bpost)
yrep <- sapply(lambda_draw, function(lam) rpois(n_sun, lam))

ppc_stats <- data.frame(
  statistic = c("mean", "variance", "maximum", "lag-1 autocorrelation"),
  observed = c(mean(y_sun_int), var(y_sun_int), max(y_sun_int), acf(y_sun_int, plot = FALSE)$acf[2]),
  ppc_median = c(median(colMeans(yrep)), median(apply(yrep, 2, var)),
                 median(apply(yrep, 2, max)), median(apply(yrep, 2, function(z) acf(z, plot = FALSE)$acf[2]))),
  ppc_lower_95 = c(quantile(colMeans(yrep), 0.025), quantile(apply(yrep, 2, var), 0.025),
                   quantile(apply(yrep, 2, max), 0.025), quantile(apply(yrep, 2, function(z) acf(z, plot = FALSE)$acf[2]), 0.025)),
  ppc_upper_95 = c(quantile(colMeans(yrep), 0.975), quantile(apply(yrep, 2, var), 0.975),
                   quantile(apply(yrep, 2, max), 0.975), quantile(apply(yrep, 2, function(z) acf(z, plot = FALSE)$acf[2]), 0.975))
)
knitr::kable(ppc_stats, digits = 3, caption = "Posterior predictive checks for the naive iid Poisson sunspot model.")
Posterior predictive checks for the naive iid Poisson sunspot model.
statistic observed ppc_median ppc_lower_95 ppc_upper_95
mean 48.633 48.647 47.576 49.779
variance 1559.101 48.738 40.851 57.455
maximum 190.000 69.000 65.000 77.000
lag-1 autocorrelation 0.814 -0.003 -0.114 0.110 
sun_df <- data.frame(year = year_sun, sunspots = y_sun)

ggplot(sun_df, aes(year, sunspots)) +
  geom_line(linewidth = 0.8) +
  labs(title = "Yearly sunspot numbers",
       subtitle = "The pattern is cyclic and dependent, not iid Poisson noise",
       x = "year", y = "sunspot number")

ppc_var <- data.frame(simulated_variance = apply(yrep, 2, var))

ggplot(ppc_var, aes(simulated_variance)) +
  geom_histogram(bins = 35, fill = "grey75", colour = "white") +
  geom_vline(xintercept = var(y_sun_int), linewidth = 1.2, linetype = 2) +
  labs(title = "Posterior predictive check: variance under iid Poisson model",
       subtitle = "Dashed line = observed variance. Large mismatch means model inadequacy.",
       x = "replicated-series variance", y = "frequency")

Key lesson. Bayesian inference is not magic. A posterior is conditional on a model. If the model ignores a major scientific structure such as periodicity, dependence or heterogeneity, posterior predictive checks should reveal the inadequacy.

7 Robust Bayesian regression: when Gaussian likelihood is misspecified

7.1 Astronomical toy motivation

Suppose we measure a simple Hubble-like relation:

\[ v_i = H d_i + \epsilon_i, \]

where \(d_i\) is distance and \(v_i\) is recession velocity. In practice, there may be outliers from measurement error, peculiar velocities, misclassification or calibration problems. A Gaussian likelihood can become overly sensitive to such outliers.

We compare:

  1. Gaussian likelihood with known \(\sigma\);
  2. Student-\(t\) likelihood with known \(\sigma\) and degrees of freedom \(\nu\), fitted by Metropolis–Hastings.
set.seed(2501)
n <- 45
d <- sort(runif(n, 5, 120))
H_true <- 70
sigma <- 420
v <- H_true * d + rnorm(n, 0, sigma)
out_id <- sample(seq_len(n), 5)
v[out_id] <- v[out_id] + rnorm(5, 2800, 900)
astro_reg <- data.frame(distance = d, velocity = v, outlier = seq_len(n) %in% out_id)

ggplot(astro_reg, aes(distance, velocity, colour = outlier)) +
  geom_point(size = 2.5) +
  scale_colour_manual(values = c("FALSE" = "black", "TRUE" = "red")) +
  geom_abline(intercept = 0, slope = H_true, linetype = 2) +
  labs(title = "Toy Hubble-like data with contaminated observations",
       subtitle = "Dashed line = true slope used in simulation",
       x = "distance", y = "velocity", colour = "contaminated")

7.1.1 Gaussian posterior for the slope

Assume \(v_i\mid H\sim N(Hd_i,\sigma^2)\) and \(H\sim N(m_0,V_0)\). Then

\[ V_n=\left(V_0^{-1}+\sigma^{-2}\sum_i d_i^2\right)^{-1}, \qquad m_n=V_n\left(V_0^{-1}m_0+\sigma^{-2}\sum_i d_i v_i\right). \]

m0 <- 70
V0 <- 30^2
Vn <- 1 / (1/V0 + sum(d^2)/sigma^2)
mn <- Vn * (m0/V0 + sum(d * v)/sigma^2)
H_gauss <- rnorm(20000, mn, sqrt(Vn))
ci_gauss <- quantile(H_gauss, c(0.025, 0.5, 0.975))

7.1.2 Student-t robust posterior by Metropolis–Hastings

The Student-\(t\) likelihood is

\[ p(v_i\mid H) = \frac{1}{\sigma}t_\nu\left(\frac{v_i-Hd_i}{\sigma}\right), \]

where \(t_\nu\) is the standard Student-\(t\) density. Heavy tails reduce the impact of a few extreme points.

logpost_t <- function(H, nu = 4) {
  sum(dt((v - H * d) / sigma, df = nu, log = TRUE) - log(sigma)) +
    dnorm(H, mean = m0, sd = sqrt(V0), log = TRUE)
}

rw_mh <- function(niter = 30000, init = 60, prop_sd = 1.2) {
  out <- numeric(niter)
  out[1] <- init
  lp <- logpost_t(init)
  acc <- 0
  for (i in 2:niter) {
    cand <- rnorm(1, out[i - 1], prop_sd)
    lp_cand <- logpost_t(cand)
    if (log(runif(1)) < lp_cand - lp) {
      out[i] <- cand
      lp <- lp_cand
      acc <- acc + 1
    } else {
      out[i] <- out[i - 1]
    }
  }
  list(draws = out, accept_rate = acc / (niter - 1))
}

mh <- rw_mh()
H_t <- mh$draws[5001:length(mh$draws)]
ci_t <- quantile(H_t, c(0.025, 0.5, 0.975))

robust_summary <- data.frame(
  model = c("Gaussian likelihood", "Student-t likelihood"),
  posterior_median = c(ci_gauss[2], ci_t[2]),
  lower_95 = c(ci_gauss[1], ci_t[1]),
  upper_95 = c(ci_gauss[3], ci_t[3]),
  true_H = H_true,
  mh_acceptance = c(NA, mh$accept_rate)
)
knitr::kable(robust_summary, digits = 2, caption = "Posterior comparison under Gaussian and robust Student-t likelihoods.")
Posterior comparison under Gaussian and robust Student-t likelihoods.
model posterior_median lower_95 upper_95 true_H mh_acceptance
Gaussian likelihood 71.62 69.94 73.31 70 NA
Student-t likelihood 69.63 67.84 71.42 70 0.63 
post_slope <- rbind(
  data.frame(H = H_gauss, model = "Gaussian"),
  data.frame(H = H_t, model = "Student-t robust")
)

ggplot(post_slope, aes(H, fill = model)) +
  geom_density(alpha = 0.45) +
  geom_vline(xintercept = H_true, linetype = 2) +
  labs(title = "Posterior for Hubble-like slope under two likelihoods",
       subtitle = "Heavy-tailed likelihood protects inference from outliers",
       x = "slope H", y = "posterior density", fill = "model")

trace_df <- data.frame(iteration = seq_along(mh$draws), H = mh$draws)

ggplot(trace_df, aes(iteration, H)) +
  geom_line(linewidth = 0.35) +
  labs(title = "Metropolis--Hastings trace for robust Student-t model",
       subtitle = paste0("Acceptance rate = ", round(mh$accept_rate, 3)),
       x = "iteration", y = "H")

Research lesson. Sometimes the likelihood is the main modelling choice. A Bayesian analysis with a badly misspecified Gaussian likelihood can still be fragile. Robust likelihoods, model checking and sensitivity analysis are part of serious Bayesian workflow.

8 Hierarchical shrinkage: borrowing strength across scientific groups

8.1 Motivation

In a basic-science setting, one may estimate a flux, reaction rate, decay rate or material property across multiple regions, instruments or experimental conditions. Some groups have many observations; others have very few. Complete pooling ignores differences; no pooling overfits small groups. Hierarchical modelling gives partial pooling.

We simulate \(J\) sky regions with true region-specific means:

\[ \mu_j\sim N(\mu_0,\tau^2), \qquad y_{ij}\mid \mu_j\sim N(\mu_j,\sigma^2). \]

set.seed(2718)
J <- 12
mu0_true <- 12
tau_true <- 3.2
sigma_y <- 2.5
n_j <- sample(2:18, J, replace = TRUE)
mu_true <- rnorm(J, mu0_true, tau_true)

y_list <- lapply(seq_len(J), function(j) rnorm(n_j[j], mu_true[j], sigma_y))
ybar <- vapply(y_list, mean, numeric(1))
se2 <- sigma_y^2 / n_j

# Empirical-Bayes estimates for hyperparameters
mu0_hat <- mean(ybar)
tau2_hat <- max(0.001, var(ybar) - mean(se2))

# Normal-normal posterior mean for each group mean
B_j <- tau2_hat / (tau2_hat + se2)
mu_shrink <- B_j * ybar + (1 - B_j) * mu0_hat

rmse_no_pool <- sqrt(mean((ybar - mu_true)^2))
rmse_shrink <- sqrt(mean((mu_shrink - mu_true)^2))

shrink_df <- data.frame(
  region = paste0("R", seq_len(J)),
  n = n_j,
  true_mean = mu_true,
  no_pooling = ybar,
  partial_pooling = mu_shrink,
  shrinkage_weight = B_j
)
knitr::kable(shrink_df, digits = 2, caption = "Hierarchical shrinkage simulation by region.")
Hierarchical shrinkage simulation by region.
region n true_mean no_pooling partial_pooling shrinkage_weight
R1 14 11.90 12.16 12.11 0.94
R2 15 10.10 10.10 10.15 0.95
R3 7 10.12 9.89 10.03 0.90
R4 13 5.55 4.45 4.85 0.94
R5 17 13.35 12.08 12.04 0.95
R6 15 13.40 12.88 12.79 0.95
R7 4 11.28 9.88 10.10 0.83
R8 2 11.61 12.09 11.84 0.71
R9 2 16.51 15.78 14.46 0.71
R10 13 12.98 13.37 13.25 0.94
R11 16 12.72 13.61 13.50 0.95
R12 5 8.71 8.23 8.65 0.86 
rmse_tab <- data.frame(
  method = c("No pooling: sample mean", "Partial pooling: hierarchical shrinkage"),
  RMSE_against_truth = c(rmse_no_pool, rmse_shrink)
)
knitr::kable(rmse_tab, digits = 3, caption = "Shrinkage usually improves unstable small-group estimation.")
Shrinkage usually improves unstable small-group estimation.
method RMSE_against_truth
No pooling: sample mean 0.769
Partial pooling: hierarchical shrinkage 0.862 
plot_df <- rbind(
  data.frame(region = seq_len(J), value = mu_true, type = "true"),
  data.frame(region = seq_len(J), value = ybar, type = "no pooling"),
  data.frame(region = seq_len(J), value = mu_shrink, type = "partial pooling")
)

ggplot(plot_df, aes(region, value, colour = type, shape = type)) +
  geom_point(size = 3) +
  geom_line(aes(group = type), linewidth = 0.6) +
  scale_x_continuous(breaks = seq_len(J)) +
  labs(title = "Hierarchical shrinkage: no pooling vs partial pooling",
       subtitle = "Partial pooling stabilizes noisy small-sample regions",
       x = "region", y = "estimated mean", colour = "quantity", shape = "quantity")

Research lesson. Hierarchical Bayesian thinking is central in natural science because measurements are often grouped: by telescope, region, detector, patient batch, material sample, or experimental run.

9 Bayes factors and the danger of over-trusting model comparison

A Bayes factor compares two models using marginal likelihoods:

\[ BF_{10}(y)=\frac{m_1(y)}{m_0(y)}. \]

If \(BF_{10}>1\), the data favour model \(M_1\) over \(M_0\). However, Bayes factors depend on the prior under each model; this is both a strength and a sensitivity issue.

9.1 Exact normal-mean Bayes factor

Suppose

\[ y_i\mid \mu\sim N(\mu,\sigma^2), \]

and we compare:

\[ M_0:\mu=0, \qquad M_1:\mu\sim N(0,\tau^2). \]

The sample mean satisfies

\[ \bar y\mid M_0\sim N(0,\sigma^2/n), \qquad \bar y\mid M_1\sim N(0,\sigma^2/n+\tau^2). \]

Therefore the Bayes factor can be computed exactly from the density of \(\bar y\) under the two models.

set.seed(3232)
n <- 30
sigma_known <- 1
mu_true <- 0.35
y_bf <- rnorm(n, mu_true, sigma_known)
ybar_bf <- mean(y_bf)

tau_grid <- seq(0.05, 3, length.out = 300)
bf10 <- dnorm(ybar_bf, 0, sqrt(sigma_known^2/n + tau_grid^2)) /
  dnorm(ybar_bf, 0, sqrt(sigma_known^2/n))
bf_df <- data.frame(tau = tau_grid, BF10 = bf10)

ggplot(bf_df, aes(tau, BF10)) +
  geom_line(linewidth = 1.1) +
  geom_hline(yintercept = 1, linetype = 2) +
  labs(title = "Bayes factor sensitivity to prior scale under the alternative",
       subtitle = paste0("Observed sample mean = ", round(ybar_bf, 3)),
       x = expression(tau), y = expression(BF[10]))

Research lesson. Bayes factors are powerful but prior-sensitive. For scientific use, report the prior scale and perform sensitivity analysis.

10 ML plus Bayesian uncertainty: galaxy velocity groups

10.1 Real data: Corona Borealis galaxy velocities

The MASS::galaxies dataset contains velocities in km/sec of 82 galaxies from six well-separated conic sections of an unfilled survey of the Corona Borealis region. Its documentation notes that multimodality in such surveys is evidence for voids and superclusters in the far universe.

Here we use ML-like mixture discovery to find possible velocity groups, and then use Bayesian normal inference to quantify the uncertainty of each group mean. This is deliberately framed as a teaching example: ML discovers structure; Bayesian inference quantifies uncertainty inside the discovered structure.

data(galaxies, package = "MASS")
gal <- as.numeric(galaxies) / 1000 # in 1000 km/s

gal_df <- data.frame(velocity = gal)

ggplot(gal_df, aes(velocity)) +
  geom_histogram(aes(y = after_stat(density)), bins = 22, fill = "grey80", colour = "white") +
  geom_density(linewidth = 1.1, colour = "blue") +
  geom_rug(alpha = 0.5) +
  labs(title = "Galaxy velocities in the Corona Borealis region",
       subtitle = "The distribution is visibly non-Gaussian and likely multimodal",
       x = "velocity (1000 km/sec)", y = "density")

10.2 Simple EM algorithm for Gaussian mixtures

We implement a small univariate Gaussian-mixture EM routine. This is not a full Bayesian mixture model, but it gives a useful bridge: it discovers latent groups; then we put Bayesian uncertainty summaries on group means.

em_mix1d <- function(x, K, nstart = 15, maxit = 500, tol = 1e-8) {
  n <- length(x)
  best <- NULL
  best_ll <- -Inf
  for (s in seq_len(nstart)) {
    if (K == 1) {
      mu <- mean(x)
      sig <- sd(x)
      pi_k <- 1
    } else {
      km <- stats::kmeans(matrix(x, ncol = 1), centers = K, nstart = 5)
      mu <- as.numeric(km$centers)
      sig <- rep(sd(x), K)
      pi_k <- as.numeric(table(factor(km$cluster, levels = 1:K))) / n
    }
    old_ll <- -Inf
    for (it in seq_len(maxit)) {
      dens <- sapply(seq_len(K), function(k) pi_k[k] * dnorm(x, mu[k], max(sig[k], 1e-4)))
      if (K == 1) dens <- matrix(dens, ncol = 1)
      row_sum <- rowSums(dens)
      row_sum[row_sum <= 0 | !is.finite(row_sum)] <- .Machine$double.eps
      resp <- dens / row_sum
      Nk <- colSums(resp)
      pi_k <- Nk / n
      mu <- colSums(resp * x) / Nk
      sig <- sqrt(colSums(resp * (x - rep(mu, each = n))^2) / Nk)
      sig <- pmax(sig, 1e-4)
      ll <- sum(log(row_sum))
      if (abs(ll - old_ll) < tol) break
      old_ll <- ll
    }
    if (ll > best_ll) {
      best_ll <- ll
      best <- list(K = K, pi = pi_k, mu = mu, sigma = sig, loglik = ll, resp = resp)
    }
  }
  p <- (K - 1) + K + K
  best$bic <- -2 * best$loglik + p * log(n)
  best
}

fits <- lapply(1:5, function(K) em_mix1d(gal, K))
bic_tab <- data.frame(
  K = 1:5,
  logLik = vapply(fits, function(z) z$loglik, numeric(1)),
  BIC = vapply(fits, function(z) z$bic, numeric(1))
)
knitr::kable(bic_tab, digits = 2, caption = "Gaussian mixture fits to galaxy velocities: smaller BIC is preferred.")
Gaussian mixture fits to galaxy velocities: smaller BIC is preferred.
K logLik BIC
1 -240.34 489.49
2 -220.06 462.15
3 -203.18 441.61
4 -202.16 452.80
5 -195.97 453.63 
best_K <- bic_tab$K[which.min(bic_tab$BIC)]
best_fit <- fits[[best_K]]
grid_g <- seq(min(gal) - 1, max(gal) + 1, length.out = 800)
comp_df <- do.call(rbind, lapply(seq_len(best_K), function(k) {
  data.frame(velocity = grid_g,
             density = best_fit$pi[k] * dnorm(grid_g, best_fit$mu[k], best_fit$sigma[k]),
             component = paste0("component ", k))
}))
assign <- max.col(best_fit$resp)
gal_df$component <- factor(assign)

ggplot() +
  geom_histogram(data = gal_df, aes(velocity, y = after_stat(density)),
                 bins = 22, fill = "grey85", colour = "white") +
  geom_line(data = comp_df, aes(velocity, density, colour = component), linewidth = 1.1) +
  geom_rug(data = gal_df, aes(velocity, colour = component), alpha = 0.7) +
  labs(title = paste("Gaussian mixture discovery for galaxy velocities: K =", best_K),
       subtitle = "ML-style latent-group discovery before Bayesian uncertainty summaries",
       x = "velocity (1000 km/sec)", y = "density", colour = "group")

10.3 Bayesian uncertainty for discovered group means

Conditional on the discovered group labels, we now compute a Bayesian normal-inverse-gamma posterior for each group mean. This is only an approximate two-stage analysis because cluster-label uncertainty is ignored. The teaching point remains valuable: after ML discovers structure, Bayesian inference can attach uncertainty to interpretable quantities.

nig_posterior <- function(y, m0 = mean(gal), k0 = 0.01, a0 = 2, b0 = 2, S = 8000) {
  n <- length(y)
  ybar <- mean(y)
  ss <- sum((y - ybar)^2)
  kn <- k0 + n
  mn <- (k0 * m0 + n * ybar) / kn
  an <- a0 + n/2
  bn <- b0 + 0.5 * ss + (k0 * n * (ybar - m0)^2) / (2 * kn)
  sigma2 <- 1 / rgamma(S, shape = an, rate = bn)
  mu <- rnorm(S, mean = mn, sd = sqrt(sigma2 / kn))
  data.frame(mu = mu, sigma = sqrt(sigma2))
}

ci_list <- lapply(levels(gal_df$component), function(cl) {
  y <- gal_df$velocity[gal_df$component == cl]
  draws <- nig_posterior(y)
  qs <- quantile(draws$mu, c(0.025, 0.5, 0.975))
  data.frame(component = cl, n = length(y), mean_observed = mean(y),
             post_median_mu = qs[2], lower_95 = qs[1], upper_95 = qs[3])
})
ci_groups <- do.call(rbind, ci_list)
knitr::kable(ci_groups, digits = 3, caption = "Bayesian uncertainty for velocity-group means after mixture discovery.")
Bayesian uncertainty for velocity-group means after mixture discovery.
component n mean_observed post_median_mu lower_95 upper_95
50% 1 3 33.044 32.999 31.492 34.439
50%1 2 72 21.400 21.401 20.896 21.912
50%2 3 7 9.710 9.729 9.094 10.369 
ci_groups$component <- factor(ci_groups$component)
ggplot(ci_groups, aes(component, post_median_mu)) +
  geom_point(size = 3) +
  geom_errorbar(aes(ymin = lower_95, ymax = upper_95), width = 0.15, linewidth = 0.9) +
  labs(title = "Posterior intervals for discovered galaxy-velocity group means",
       x = "discovered group", y = "velocity-group mean (1000 km/sec)")

Research lesson. ML is excellent for discovering patterns such as clusters. Bayesian inference is excellent for quantifying uncertainty in scientifically meaningful quantities after or during pattern discovery. A fully Bayesian mixture model would combine both steps in a single posterior, but the two-stage workflow is already useful for teaching.

11 Prior sensitivity: a necessary scientific habit

Prior sensitivity is not a weakness; it is a diagnostic. If conclusions change drastically under reasonable priors, the data may not be strong enough to settle the question.

We revisit the Beta–Binomial experiment with multiple priors.

priors <- data.frame(a = c(1, 2, 10, 30), b = c(1, 8, 10, 70), label = c("Uniform Beta(1,1)", "Low-detection Beta(2,8)", "Centered Beta(10,10)", "Strong low Beta(30,70)"))
post_sens <- do.call(rbind, lapply(seq_len(nrow(priors)), function(i) {
  a <- priors$a[i]; b <- priors$b[i]
  data.frame(theta = th,
             density = dbeta(th, a + y, b + n - y),
             prior = priors$label[i])
}))

ggplot(post_sens, aes(theta, density, colour = prior)) +
  geom_line(linewidth = 1.05) +
  geom_vline(xintercept = y/n, linetype = 2) +
  labs(title = "Prior sensitivity of the Beta--Binomial posterior",
       subtitle = "Dashed line = observed detection fraction",
       x = expression(theta), y = "posterior density", colour = "prior")

Lesson. With enough information, reasonable priors often lead to similar posteriors. With limited information, the prior is visible. That visibility is not a flaw; it is scientific honesty.

12 Compact Bayesian workflow for scientific research

workflow <- data.frame(
  Step = c("1. Scientific question", "2. Data-generating story", "3. Prior construction", "4. Posterior computation", "5. Posterior summaries", "6. Predictive checks", "7. Sensitivity", "8. Scientific decision"),
  Question = c("What is the estimand or latent object?", "What is observed and how is it noisy?", "What information or regularization is defensible?", "Can we use conjugacy, quadrature, MCMC or SMC?", "What are credible intervals, probabilities and predictions?", "Can the model reproduce key features of data?", "Do conclusions depend on priors/likelihood choices?", "What action or research conclusion follows?"),
  Example = c("source intensity, galaxy group, solar-cycle parameter", "Poisson photons, background, telescope noise", "weakly informative gamma/normal/t priors", "grid posterior, MCMC, EM + Bayes", "posterior mean, P(effect > 0), posterior predictive", "sunspot variance and autocorrelation", "Gaussian vs Student-t likelihood", "follow-up observation, classify, publish uncertainty")
)
knitr::kable(workflow, caption = "A practical Bayesian workflow for basic-science researchers.")
A practical Bayesian workflow for basic-science researchers.
Step Question Example
1. Scientific question What is the estimand or latent object? source intensity, galaxy group, solar-cycle parameter
2. Data-generating story What is observed and how is it noisy? Poisson photons, background, telescope noise
3. Prior construction What information or regularization is defensible? weakly informative gamma/normal/t priors
4. Posterior computation Can we use conjugacy, quadrature, MCMC or SMC? grid posterior, MCMC, EM + Bayes
5. Posterior summaries What are credible intervals, probabilities and predictions? posterior mean, P(effect > 0), posterior predictive
6. Predictive checks Can the model reproduce key features of data? sunspot variance and autocorrelation
7. Sensitivity Do conclusions depend on priors/likelihood choices? Gaussian vs Student-t likelihood
8. Scientific decision What action or research conclusion follows? follow-up observation, classify, publish uncertainty

13 Final message

Bayesian inference is not merely a formula. It is a probability environment for scientific learning:

\[ \text{model} + \text{prior information} + \text{data} \longrightarrow \text{posterior uncertainty} + \text{predictive checking} + \text{decision}. \]

The examples above show several important principles:

14 References and further reading