1 Aim of this note

This note is a standalone RPubs lecture note on Bayesian statistical inference for students and researchers in the basic sciences. The emphasis is not on black-box Bayesian software. The aim is to understand the probability environment, prior, likelihood, posterior, posterior prediction, uncertainty quantification, simulation verification, and model checking.

The scientific flavor is deliberately astronomical:

  • a toy Hubble-law example for estimating an expansion-like slope under noisy velocity measurements;
  • solar activity using the built-in sunspot.year time series;
  • galaxy velocities using the MASS::galaxies dataset;
  • simulation examples that show how Bayesian inference behaves when the model is correct and how model checking warns us when the model is wrong.

Core message. Bayesian inference is probabilistic inference: unknown scientific quantities are represented through a probability distribution, updated by data, and then used for estimation, prediction, and decision-making under uncertainty.

2 Probability environment: what is random, what is observed?

A probability model begins with a probability space

\[ (\Omega,\mathcal F,P), \]

where:

  • \(\Omega\) is the sample space of elementary outcomes;
  • \(\mathcal F\) is a sigma-field of events;
  • \(P\) is a probability measure satisfying \(P(\Omega)=1\).

A random variable is a measurable function

\[ Y:\Omega\to \mathcal Y. \]

In statistics, we observe a realized value \(y\) of \(Y\), and then infer something about an unknown parameter \(\theta\) or latent scientific quantity.

A statistical model is a family of probability laws

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

where \(\Theta\) is the parameter space. In an astronomical example, \(\theta\) could be a solar-activity rate, a galaxy-cluster velocity center, an expansion slope, or a latent state in a filtering problem.

2.1 Frequentist and Bayesian environments

In a frequentist statistical environment, \(\theta\) is fixed but unknown. Probability statements are made about repeated samples \(Y\sim P_\theta\).

In a Bayesian statistical environment, \(\theta\) is also uncertain and is assigned a prior probability measure \(\Pi\) on \((\Theta,\mathcal T)\):

\[ \theta\sim \Pi. \]

Then data are generated conditionally:

\[ Y\mid \theta \sim P_\theta. \]

After observing \(Y=y\), we update from prior to posterior.

3 Likelihood, prior, posterior

Suppose \(Y\) has density or probability mass function \(p_\theta(y)\). The likelihood is

\[ L(\theta;y)=p_\theta(y), \]

viewed as a function of \(\theta\) after observing \(y\). The likelihood is not a probability distribution over \(\theta\) by itself.

If the prior has density \(\pi(\theta)\), Bayes’ formula gives

\[ \pi(\theta\mid y) = \frac{p_\theta(y)\pi(\theta)}{\int_\Theta p_u(y)\pi(u)\,du}. \]

The denominator

\[ m(y)=\int_\Theta p_\theta(y)\pi(\theta)\,d\theta \]

is called the marginal likelihood or evidence.

For a measurable set \(B\subseteq \Theta\), the posterior probability is

\[ \Pi(B\mid Y=y) = \frac{\int_B p_\theta(y)\,\Pi(d\theta)} {\int_\Theta p_\theta(y)\,\Pi(d\theta)}. \]

Interpretation. A posterior statement such as \(P(\theta>0\mid y)=0.98\) means: under the assumed model and prior, after observing the data, 98% of the posterior probability mass lies on positive values of \(\theta\).

4 Posterior summaries

After obtaining \(\pi(\theta\mid y)\), common Bayesian summaries are:

Quantity Definition Meaning
Posterior mean \(E(\theta\mid y)\) squared-error optimal estimate
Posterior median median of \(\pi(\theta\mid y)\) robust central estimate
MAP \(\arg\max_\theta \pi(\theta\mid y)\) posterior mode
Credible interval interval \(C(y)\) with \(\Pi(C(y)\mid y)=0.95\) uncertainty interval for \(\theta\)
Posterior tail probability \(P(\theta>c\mid y)\) probability of scientific statement
Posterior predictive \(p(\tilde y\mid y)\) uncertainty for future or replicated data

The posterior predictive distribution is

\[ p(\tilde y\mid y)=\int p(\tilde y\mid \theta)\pi(\theta\mid y)\,d\theta. \]

This is central in scientific work because science often needs prediction with uncertainty, not only parameter estimation.

5 Simulation 1: Beta–Binomial updating

Suppose \(Y\mid p\sim\text{Binomial}(n,p)\), where \(p\) is an unknown probability. Use the prior

\[ p\sim\text{Beta}(a,b). \]

If \(Y=y\), the posterior is

\[ p\mid y\sim \text{Beta}(a+y,b+n-y). \]

This is a conjugate model: prior and posterior belong to the same family.

# True data-generating probability
p_true <- 0.63
n <- 40
y <- rbinom(1, size = n, prob = p_true)

# Prior Beta(a,b)
a <- 2
b <- 2

# Posterior Beta(a+y, b+n-y)
a_post <- a + y
b_post <- b + n - y

p_grid <- seq(0.001, 0.999, length.out = 1000)
prior <- dbeta(p_grid, a, b)
lik <- dbinom(y, size = n, prob = p_grid)
lik_scaled <- lik / max(lik) * max(prior)
post <- dbeta(p_grid, a_post, b_post)

df <- data.frame(
  p = rep(p_grid, 3),
  density = c(prior, lik_scaled, post),
  curve = rep(c("Prior Beta(2,2)", "Scaled likelihood", "Posterior"), each = length(p_grid))
)

ggplot(df, aes(p, density, color = curve)) +
  geom_line(linewidth = 1.2) +
  geom_vline(xintercept = p_true, linetype = "dashed", color = "black") +
  labs(
    title = "Bayesian updating in a Binomial experiment",
    subtitle = paste0("Observed y = ", y, " successes out of n = ", n, "; dashed line is true p = ", p_true),
    x = "p", y = "density / scaled likelihood", color = ""
  )

Posterior numerical summary:

beta_summary <- data.frame(
  quantity = c("observed successes", "sample proportion", "posterior mean", "posterior median", "95% credible lower", "95% credible upper", "P(p > 0.5 | data)"),
  value = c(
    y,
    y/n,
    a_post/(a_post+b_post),
    qbeta(0.5, a_post, b_post),
    qbeta(0.025, a_post, b_post),
    qbeta(0.975, a_post, b_post),
    1 - pbeta(0.5, a_post, b_post)
  )
)
knitr::kable(beta_summary, digits = 4)
quantity value
observed successes 28.0000
sample proportion 0.7000
posterior mean 0.6818
posterior median 0.6846
95% credible lower 0.5387
95% credible upper 0.8092
P(p > 0.5 | data) 0.9931

5.1 Verification by repeated simulation

A Bayesian credible interval is not defined by repeated-sampling coverage, but repeated simulation is still useful for understanding behavior. Here we simulate many datasets from the same true \(p\), compute the posterior 95% interval, and check how often it contains the true value.

set.seed(12)
B <- 1000
n <- 40
p_true <- 0.63
cover <- numeric(B)
post_mean <- numeric(B)

for (r in seq_len(B)) {
  yr <- rbinom(1, n, p_true)
  lo <- qbeta(0.025, a + yr, b + n - yr)
  hi <- qbeta(0.975, a + yr, b + n - yr)
  cover[r] <- (lo <= p_true && p_true <= hi)
  post_mean[r] <- (a + yr) / (a + b + n)
}

res <- data.frame(
  simulation_quantity = c("empirical containment of true p", "mean posterior mean", "SD of posterior means"),
  value = c(mean(cover), mean(post_mean), sd(post_mean))
)
knitr::kable(res, digits = 4)
simulation_quantity value
empirical containment of true p 0.9570
mean posterior mean 0.6196
SD of posterior means 0.0675 
hist_df <- data.frame(post_mean = post_mean)
ggplot(hist_df, aes(post_mean)) +
  geom_histogram(bins = 35, fill = "steelblue", color = "white") +
  geom_vline(xintercept = p_true, color = "red", linewidth = 1.1) +
  labs(
    title = "Repeated simulation: posterior mean across experiments",
    subtitle = "Red line is the true p. The posterior mean varies because the data vary.",
    x = "posterior mean", y = "count"
  )

6 Simulation 2: Normal mean with known variance

Suppose

\[ Y_i\mid \mu \sim N(\mu,\sigma^2),\qquad i=1,\ldots,n, \]

where \(\sigma\) is known. Use prior

\[ \mu\sim N(\mu_0,\tau_0^2). \]

Then

\[ \mu\mid y\sim N(\mu_n,\tau_n^2), \]

where

\[ \tau_n^2=\left(\frac{1}{\tau_0^2}+\frac{n}{\sigma^2}\right)^{-1}, \qquad \mu_n=\tau_n^2\left(\frac{\mu_0}{\tau_0^2}+\frac{n\bar y}{\sigma^2}\right). \]

This formula shows shrinkage: the posterior mean is a precision-weighted compromise between prior mean and sample mean.

set.seed(7)
mu_true <- 2.4
sigma <- 1.2
n <- 25
y <- rnorm(n, mu_true, sigma)

mu0 <- 0
tau0 <- 3

tau_n2 <- 1 / (1/tau0^2 + n/sigma^2)
mu_n <- tau_n2 * (mu0/tau0^2 + n*mean(y)/sigma^2)

mu_grid <- seq(-1.5, 4.5, length.out = 1000)
prior <- dnorm(mu_grid, mu0, tau0)
post <- dnorm(mu_grid, mu_n, sqrt(tau_n2))
lik_scaled <- dnorm(mu_grid, mean(y), sigma/sqrt(n))
lik_scaled <- lik_scaled / max(lik_scaled) * max(post)

ndf <- data.frame(
  mu = rep(mu_grid, 3),
  density = c(prior, lik_scaled, post),
  curve = rep(c("Prior", "Sampling information (scaled)", "Posterior"), each = length(mu_grid))
)

ggplot(ndf, aes(mu, density, color = curve)) +
  geom_line(linewidth = 1.2) +
  geom_vline(xintercept = mu_true, linetype = "dashed", color = "black") +
  labs(
    title = "Normal mean inference: prior + likelihood = posterior",
    subtitle = "The posterior shrinks the sample mean toward the prior mean according to precision.",
    x = expression(mu), y = "density", color = ""
  )

summary_normal <- data.frame(
  quantity = c("sample mean", "prior mean", "posterior mean", "posterior SD", "95% lower", "95% upper"),
  value = c(mean(y), mu0, mu_n, sqrt(tau_n2), qnorm(0.025, mu_n, sqrt(tau_n2)), qnorm(0.975, mu_n, sqrt(tau_n2)))
)
knitr::kable(summary_normal, digits = 4)
quantity value
sample mean 2.9568
prior mean 0.0000
posterior mean 2.9380
posterior SD 0.2392
95% lower 2.4691
95% upper 3.4069

7 Simulation 3: Astronomical toy example — estimating a Hubble-like slope

A very simple astronomical model says that recession velocity \(v_i\) is approximately proportional to distance \(d_i\):

\[ v_i = H d_i + \epsilon_i, \qquad \epsilon_i\sim N(0,\sigma^2). \]

Here \(H\) is a slope parameter. In cosmology this resembles Hubble’s law, but here it is only a toy demonstration.

Let

\[ H\sim N(H_0,s_0^2). \]

With known \(\sigma\), the posterior is again normal:

\[ H\mid v,d \sim N(m_n,s_n^2), \]

where

\[ s_n^2=\left(\frac{1}{s_0^2}+\frac{\sum_i d_i^2}{\sigma^2}\right)^{-1}, \qquad m_n=s_n^2\left(\frac{H_0}{s_0^2}+\frac{\sum_i d_i v_i}{\sigma^2}\right). \]

set.seed(2026)
n <- 45
H_true <- 70
sigma_v <- 250
D <- sort(runif(n, 5, 180))
V <- H_true * D + rnorm(n, 0, sigma_v)

H0 <- 65
s0 <- 20
sn2 <- 1 / (1/s0^2 + sum(D^2)/sigma_v^2)
mn <- sn2 * (H0/s0^2 + sum(D*V)/sigma_v^2)

ols_slope <- sum(D*V) / sum(D^2)

H_grid <- seq(45, 95, length.out = 1000)
prior_H <- dnorm(H_grid, H0, s0)
post_H <- dnorm(H_grid, mn, sqrt(sn2))

dh <- data.frame(distance = D, velocity = V)

p1 <- ggplot(dh, aes(distance, velocity)) +
  geom_point(size = 2.3, alpha = 0.8) +
  geom_abline(intercept = 0, slope = H_true, color = "black", linetype = "dashed", linewidth = 1) +
  geom_abline(intercept = 0, slope = ols_slope, color = "orange", linewidth = 1) +
  geom_abline(intercept = 0, slope = mn, color = "blue", linewidth = 1) +
  labs(
    title = "Toy astronomical velocity--distance relation",
    subtitle = "Black dashed: true slope; orange: least squares; blue: Bayesian posterior mean.",
    x = "distance (toy units)", y = "recession velocity (toy units)"
  )
print(p1)

p2 <- ggplot(data.frame(H = rep(H_grid, 2), density = c(prior_H, post_H), curve = rep(c("prior", "posterior"), each = length(H_grid))),
             aes(H, density, color = curve)) +
  geom_line(linewidth = 1.2) +
  geom_vline(xintercept = H_true, linetype = "dashed") +
  labs(title = "Posterior uncertainty for the expansion-like slope", x = "H", y = "density", color = "")
print(p2)

h_summary <- data.frame(
  quantity = c("true H", "least-squares slope", "posterior mean", "posterior SD", "95% lower", "95% upper"),
  value = c(H_true, ols_slope, mn, sqrt(sn2), qnorm(0.025, mn, sqrt(sn2)), qnorm(0.975, mn, sqrt(sn2)))
)
knitr::kable(h_summary, digits = 3)
quantity value
true H 70.000
least-squares slope 69.789
posterior mean 69.787
posterior SD 0.437
95% lower 68.931
95% upper 70.643

8 Model misspecification: outliers and robust Bayesian likelihoods

Scientific measurements can be contaminated by peculiar velocities, calibration problems, or wrong identification. If the likelihood assumes perfect Gaussian noise, an outlier can strongly affect the posterior.

We compare two likelihoods for the same slope \(H\):

\[ v_i\mid H \sim N(Hd_i,\sigma^2), \]

and a heavier-tailed Student likelihood

\[ \frac{v_i-Hd_i}{\sigma}\sim t_\nu. \]

The Student likelihood is less surprised by extreme residuals and therefore more robust.

set.seed(2027)
V_bad <- V
V_bad[which.max(D)] <- V_bad[which.max(D)] + 2600

H_grid <- seq(40, 105, length.out = 1200)
log_prior <- dnorm(H_grid, H0, s0, log = TRUE)

loglik_normal <- sapply(H_grid, function(H) sum(dnorm(V_bad, mean = H * D, sd = sigma_v, log = TRUE)))
loglik_t <- sapply(H_grid, function(H) sum(dt((V_bad - H * D)/sigma_v, df = 3, log = TRUE) - log(sigma_v)))

normalize_logpost <- function(logp) {
  w <- exp(logp - max(logp))
  w / sum(w)
}

post_norm_w <- normalize_logpost(log_prior + loglik_normal)
post_t_w <- normalize_logpost(log_prior + loglik_t)

grid_mean <- function(x, w) sum(x * w)
grid_quantile <- function(x, w, probs) {
  cw <- cumsum(w)
  sapply(probs, function(p) x[which(cw >= p)[1]])
}

mn_norm <- grid_mean(H_grid, post_norm_w)
mn_t <- grid_mean(H_grid, post_t_w)
ci_norm <- grid_quantile(H_grid, post_norm_w, c(0.025, 0.975))
ci_t <- grid_quantile(H_grid, post_t_w, c(0.025, 0.975))

p_bad <- ggplot(data.frame(distance = D, velocity = V_bad), aes(distance, velocity)) +
  geom_point(size = 2.4, alpha = 0.8) +
  geom_abline(intercept = 0, slope = H_true, linetype = "dashed", color = "black", linewidth = 1) +
  geom_abline(intercept = 0, slope = mn_norm, color = "red", linewidth = 1) +
  geom_abline(intercept = 0, slope = mn_t, color = "blue", linewidth = 1) +
  labs(
    title = "Outlier-contaminated velocity--distance data",
    subtitle = "Red: Gaussian posterior mean; blue: robust Student-t posterior mean; dashed: true slope.",
    x = "distance", y = "velocity"
  )
print(p_bad)

post_df <- data.frame(
  H = rep(H_grid, 2),
  posterior_mass = c(post_norm_w, post_t_w),
  model = rep(c("Gaussian likelihood", "Student-t likelihood"), each = length(H_grid))
)

ggplot(post_df, aes(H, posterior_mass, color = model)) +
  geom_line(linewidth = 1.2) +
  geom_vline(xintercept = H_true, linetype = "dashed") +
  labs(
    title = "Posterior under model misspecification",
    subtitle = "Heavy-tailed likelihood protects inference from one extreme observation.",
    x = "H", y = "posterior grid mass", color = ""
  )

robust_table <- data.frame(
  model = c("Gaussian likelihood", "Student-t likelihood"),
  posterior_mean = c(mn_norm, mn_t),
  lower_95 = c(ci_norm[1], ci_t[1]),
  upper_95 = c(ci_norm[2], ci_t[2])
)
knitr::kable(robust_table, digits = 3)
model posterior_mean lower_95 upper_95
Gaussian likelihood 71.091 70.250 71.931
Student-t likelihood 69.997 68.949 71.009

Lesson. Bayesian inference is not automatically safe. It is conditional on the likelihood and prior. A wrong likelihood can produce confidently wrong conclusions. Robust Bayes changes the probability model so that unusual observations do not dominate the inference.

9 Real astronomical data 1: yearly sunspot numbers

The built-in R dataset sunspot.year contains yearly sunspot numbers from 1700 to 1988 and has 289 observations. Sunspots are a classical observable of solar activity.

sun <- data.frame(
  year = as.numeric(time(sunspot.year)),
  count = as.numeric(sunspot.year)
)

knitr::kable(head(sun, 8), caption = "First few observations of yearly sunspot numbers")
First few observations of yearly sunspot numbers
year count
1700 5
1701 11
1702 16
1703 23
1704 36
1705 58
1706 29
1707 20 
ggplot(sun, aes(year, count)) +
  geom_line(color = "darkorange", linewidth = 0.8) +
  labs(
    title = "Yearly sunspot numbers, 1700--1988",
    subtitle = "Solar activity is cyclic and nonstationary; a simple iid model is only a first approximation.",
    x = "year", y = "sunspot number"
  )

9.1 A simple Gamma–Poisson rate comparison

For a count \(Y\), a basic model is

\[ Y_i\mid \lambda \sim \text{Poisson}(\lambda), \qquad \lambda\sim\text{Gamma}(a,b), \]

where the Gamma distribution is parameterized by shape \(a\) and rate \(b\). If \(y_1,\ldots,y_n\) are observed, then

\[ \lambda\mid y\sim \text{Gamma}\left(a+\sum_i y_i,\, b+n\right). \]

We compare average activity before and after 1850. This is deliberately simple and should not be mistaken for a serious solar-cycle model.

early <- subset(sun, year < 1850)$count
late <- subset(sun, year >= 1850)$count

a0 <- 1
b0 <- 0.05

post_early <- c(shape = a0 + sum(early), rate = b0 + length(early))
post_late <- c(shape = a0 + sum(late), rate = b0 + length(late))

S <- 20000
lambda_early <- rgamma(S, shape = post_early["shape"], rate = post_early["rate"])
lambda_late <- rgamma(S, shape = post_late["shape"], rate = post_late["rate"])

sun_table <- data.frame(
  period = c("1700--1849", "1850--1988"),
  n_years = c(length(early), length(late)),
  sample_mean = c(mean(early), mean(late)),
  posterior_mean_rate = c(mean(lambda_early), mean(lambda_late)),
  lower_95 = c(quantile(lambda_early, 0.025), quantile(lambda_late, 0.025)),
  upper_95 = c(quantile(lambda_early, 0.975), quantile(lambda_late, 0.975))
)
knitr::kable(sun_table, digits = 3)
period n_years sample_mean posterior_mean_rate lower_95 upper_95
1700–1849 150 43.795 43.781 42.720 44.864
1850–1988 139 53.814 53.799 52.574 55.031 
p_late_gt_early <- mean(lambda_late > lambda_early)
cat("Posterior probability that late-period mean sunspot rate exceeds early-period rate:", round(p_late_gt_early, 4), "\n")
## Posterior probability that late-period mean sunspot rate exceeds early-period rate: 1
rate_df <- data.frame(
  lambda = c(lambda_early, lambda_late),
  period = rep(c("1700--1849", "1850--1988"), each = S)
)

ggplot(rate_df, aes(lambda, fill = period, color = period)) +
  geom_density(alpha = 0.25, linewidth = 1.1) +
  labs(
    title = "Posterior distributions for mean yearly sunspot rate",
    subtitle = "Gamma--Poisson model: useful as a first uncertainty calculation, but too simple for solar cycles.",
    x = expression(lambda), y = "posterior density", fill = "", color = ""
  )

9.2 Posterior predictive check: why the simple Poisson model is inadequate

If \(Y_i\sim \text{Poisson}(\lambda)\), then \(\operatorname{Var}(Y_i)\approx E(Y_i)\). Sunspot data have cyclicity, dependence, and overdispersion. A posterior predictive check makes this visible.

y_all <- sun$count
n_all <- length(y_all)
post_all <- c(shape = a0 + sum(y_all), rate = b0 + n_all)
lambda_all <- rgamma(2000, shape = post_all["shape"], rate = post_all["rate"])

var_to_mean <- function(x) var(x) / mean(x)
obs_vtm <- var_to_mean(y_all)
ppc_vtm <- sapply(lambda_all, function(lam) var_to_mean(rpois(n_all, lam)))

ppc_df <- data.frame(vtm = ppc_vtm)
ggplot(ppc_df, aes(vtm)) +
  geom_histogram(bins = 40, fill = "skyblue", color = "white") +
  geom_vline(xintercept = obs_vtm, color = "red", linewidth = 1.2) +
  labs(
    title = "Posterior predictive check for Poisson sunspot model",
    subtitle = "Red line is observed variance/mean ratio. The simple iid Poisson model misses solar-cycle variability.",
    x = "variance / mean", y = "posterior predictive frequency"
  )

ppc_sun_table <- data.frame(
  statistic = c("observed variance/mean", "posterior predictive mean", "posterior predictive 95% lower", "posterior predictive 95% upper", "P(T_rep >= T_obs)"),
  value = c(obs_vtm, mean(ppc_vtm), quantile(ppc_vtm, 0.025), quantile(ppc_vtm, 0.975), mean(ppc_vtm >= obs_vtm))
)
knitr::kable(ppc_sun_table, digits = 4)
statistic value
observed variance/mean 32.0529
posterior predictive mean 1.0000
posterior predictive 95% lower 0.8482
posterior predictive 95% upper 1.1712
P(T_rep >= T_obs) 0.0000

Scientific lesson. A Bayesian posterior can be mathematically correct for the model and still scientifically inadequate if the model ignores time dependence. This is why posterior predictive checking is an essential part of Bayesian workflow.

10 Real astronomical data 2: galaxy velocities and posterior predictive checking

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

gal <- as.numeric(MASS::galaxies) / 1000  # thousand km/sec
gal_df <- data.frame(velocity = gal)

summary_gal <- data.frame(
  n = length(gal),
  mean = mean(gal),
  sd = sd(gal),
  median = median(gal),
  min = min(gal),
  max = max(gal)
)
knitr::kable(summary_gal, digits = 3, caption = "Galaxy velocities in 1000 km/sec")
Galaxy velocities in 1000 km/sec
n mean sd median min max
82 20.828 4.564 20.834 9.172 34.279 
ggplot(gal_df, aes(velocity)) +
  geom_histogram(aes(y = after_stat(density)), bins = 18, fill = "gray85", color = "white") +
  geom_density(color = "darkblue", linewidth = 1.2) +
  geom_rug(alpha = 0.6) +
  labs(
    title = "Galaxy velocities: visual evidence of non-Gaussian structure",
    subtitle = "A single normal distribution is likely too simple for clustered cosmic structure.",
    x = "velocity (1000 km/sec)", y = "density"
  )

10.1 Bayesian normal model for galaxy velocities

As a first model, suppose

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

Use the conjugate Normal–Inverse-Gamma prior:

\[ \mu\mid \sigma^2\sim N(\mu_0,\sigma^2/\kappa_0), \qquad \sigma^2\sim \text{Inv-Gamma}(\alpha_0,\beta_0). \]

nig_posterior <- function(y, mu0 = 20, kappa0 = 0.01, alpha0 = 2, beta0 = 20) {
  n <- length(y)
  ybar <- mean(y)
  ss <- sum((y - ybar)^2)
  kappa_n <- kappa0 + n
  mu_n <- (kappa0 * mu0 + n * ybar) / kappa_n
  alpha_n <- alpha0 + n / 2
  beta_n <- beta0 + 0.5 * ss + (kappa0 * n * (ybar - mu0)^2) / (2 * kappa_n)
  list(mu_n = mu_n, kappa_n = kappa_n, alpha_n = alpha_n, beta_n = beta_n)
}

sample_nig <- function(post, S = 5000) {
  sigma2 <- 1 / rgamma(S, shape = post$alpha_n, rate = post$beta_n)
  mu <- rnorm(S, mean = post$mu_n, sd = sqrt(sigma2 / post$kappa_n))
  data.frame(mu = mu, sigma = sqrt(sigma2))
}

max_gap <- function(x) max(diff(sort(x)))
post_gal <- nig_posterior(gal)
draw_gal <- sample_nig(post_gal, S = 8000)

post_sum_gal <- data.frame(
  parameter = c("mu", "sigma"),
  posterior_mean = c(mean(draw_gal$mu), mean(draw_gal$sigma)),
  lower_95 = c(quantile(draw_gal$mu, 0.025), quantile(draw_gal$sigma, 0.025)),
  upper_95 = c(quantile(draw_gal$mu, 0.975), quantile(draw_gal$sigma, 0.975))
)
knitr::kable(post_sum_gal, digits = 3)
parameter posterior_mean lower_95 upper_95
mu 20.819 19.818 21.781
sigma 4.520 3.893 5.278 
mu_df <- data.frame(mu = draw_gal$mu)
ggplot(mu_df, aes(mu)) +
  geom_density(fill = "lightgreen", alpha = 0.5, color = "darkgreen", linewidth = 1.1) +
  labs(
    title = "Posterior distribution for the mean galaxy velocity under a single-normal model",
    x = expression(mu~"(1000 km/sec)"), y = "posterior density"
  )

10.2 Posterior predictive check for cosmic clustering

A single normal model can estimate a mean and variance, but it may not reproduce multimodality. We use the largest gap between sorted velocities as a simple discrepancy statistic:

\[ T(y)=\max_i \{y_{(i+1)}-y_{(i)}\}. \]

Large gaps suggest separated clusters.

set.seed(19)
Spp <- 2000
idx <- sample(seq_len(nrow(draw_gal)), Spp)
pp_gaps <- numeric(Spp)
for (s in seq_len(Spp)) {
  yy <- rnorm(length(gal), mean = draw_gal$mu[idx[s]], sd = draw_gal$sigma[idx[s]])
  pp_gaps[s] <- max_gap(yy)
}
obs_gap <- max_gap(gal)

pp_gap_df <- data.frame(max_gap = pp_gaps)
ggplot(pp_gap_df, aes(max_gap)) +
  geom_histogram(bins = 40, fill = "lavender", color = "white") +
  geom_vline(xintercept = obs_gap, color = "red", linewidth = 1.2) +
  labs(
    title = "Posterior predictive check for a single-normal galaxy model",
    subtitle = "Red line: observed largest velocity gap. Tail behavior suggests model inadequacy if red is extreme.",
    x = "largest gap between sorted velocities", y = "posterior predictive frequency"
  )

gap_table <- data.frame(
  statistic = c("observed largest gap", "posterior predictive mean gap", "posterior predictive 95% lower", "posterior predictive 95% upper", "Bayesian p-value P(T_rep >= T_obs)"),
  value = c(obs_gap, mean(pp_gaps), quantile(pp_gaps, 0.025), quantile(pp_gaps, 0.975), mean(pp_gaps >= obs_gap))
)
knitr::kable(gap_table, digits = 4)
statistic value
observed largest gap 5.6780
posterior predictive mean gap 2.8984
posterior predictive 95% lower 1.1520
posterior predictive 95% upper 6.5742
Bayesian p-value P(T_rep >= T_obs) 0.0485

10.3 Exploratory mixture thinking: ML clustering plus Bayesian uncertainty

A full Bayesian mixture model is beyond this introductory note. But we can show a useful hybrid idea:

  1. use an exploratory algorithm to identify possible clusters;
  2. do Bayesian uncertainty quantification within each cluster.

This is not a replacement for a full Bayesian mixture, but it is pedagogically useful.

set.seed(25)
km <- kmeans(gal, centers = 3, nstart = 50)
cl_means <- tapply(gal, km$cluster, mean)
ord <- order(cl_means)
cluster_ordered <- match(km$cluster, ord)

gal_cl <- data.frame(
  velocity = gal,
  cluster = factor(cluster_ordered, labels = c("low-velocity group", "middle group", "high-velocity group"))
)

ggplot(gal_cl, aes(velocity, fill = cluster)) +
  geom_histogram(position = "identity", alpha = 0.55, bins = 18, color = "white") +
  geom_rug(aes(color = cluster), alpha = 0.7) +
  labs(
    title = "Exploratory clustering of galaxy velocities",
    subtitle = "ML-like clustering proposes groups; Bayesian inference can quantify group centers and uncertainty.",
    x = "velocity (1000 km/sec)", y = "count", fill = "", color = ""
  )

cluster_levels <- levels(gal_cl$cluster)
cluster_results <- do.call(rbind, lapply(cluster_levels, function(cl) {
  ycl <- gal_cl$velocity[gal_cl$cluster == cl]
  pst <- nig_posterior(ycl, mu0 = mean(gal), kappa0 = 0.01, alpha0 = 2, beta0 = 20)
  dr <- sample_nig(pst, S = 6000)
  data.frame(
    cluster = cl,
    n = length(ycl),
    sample_mean = mean(ycl),
    posterior_mean_center = mean(dr$mu),
    lower_95_center = quantile(dr$mu, 0.025),
    upper_95_center = quantile(dr$mu, 0.975)
  )
}))

knitr::kable(cluster_results, digits = 3, caption = "Bayesian uncertainty for exploratory galaxy velocity groups")
Bayesian uncertainty for exploratory galaxy velocity groups
cluster n sample_mean posterior_mean_center lower_95_center upper_95_center
2.5% low-velocity group 7 9.710 9.713 8.041 11.318
2.5%1 middle group 70 21.245 21.242 20.727 21.742
2.5%2 high-velocity group 5 30.564 30.538 27.300 33.811

Lesson for scientific data science. ML can help detect structure; Bayesian inference can quantify uncertainty about that structure. In scientific work, these two are often complementary rather than competitors.

11 Prior choice: subjective, objective, weakly informative

A prior is not merely a mathematical nuisance. It expresses a probability distribution on unknown quantities before the current data are used.

Prior type Meaning Scientific use
Subjective prior expresses expert knowledge or previous experiments instrument calibration, known physical ranges
Objective/reference prior designed from model structure, often to be weakly informative default inference when little prior information is available
Weakly informative prior rules out absurd values while remaining broad stabilizes high-dimensional or weak-data problems
Conjugate prior algebraically convenient prior family transparent teaching and fast computation
Hierarchical prior parameters share information through common hyperparameters partial pooling across telescopes, regions, species, experiments

11.1 Prior sensitivity demonstration

Small data can be sensitive to prior choice. Large data usually reduce this sensitivity, unless the likelihood is misspecified.

set.seed(99)
n_small <- 8
p_true <- 0.7
y_small <- rbinom(1, n_small, p_true)
priors <- data.frame(
  prior = c("Uniform Beta(1,1)", "Skeptical Beta(2,8)", "Optimistic Beta(8,2)", "Moderate Beta(3,3)"),
  a = c(1, 2, 8, 3),
  b = c(1, 8, 2, 3)
)

prior_sens <- priors
prior_sens$post_mean <- (priors$a + y_small) / (priors$a + priors$b + n_small)
prior_sens$lower_95 <- qbeta(0.025, priors$a + y_small, priors$b + n_small - y_small)
prior_sens$upper_95 <- qbeta(0.975, priors$a + y_small, priors$b + n_small - y_small)
prior_sens$prob_gt_half <- 1 - pbeta(0.5, priors$a + y_small, priors$b + n_small - y_small)

knitr::kable(prior_sens, digits = 3, caption = paste0("Prior sensitivity with small data: y = ", y_small, " out of n = ", n_small))
Prior sensitivity with small data: y = 5 out of n = 8
prior a b post_mean lower_95 upper_95 prob_gt_half
Uniform Beta(1,1) 1 1 0.600 0.299 0.863 0.746
Skeptical Beta(2,8) 2 8 0.389 0.184 0.617 0.166
Optimistic Beta(8,2) 8 2 0.722 0.501 0.897 0.975
Moderate Beta(3,3) 3 3 0.571 0.316 0.808 0.709 
plot_prior_sens <- do.call(rbind, lapply(seq_len(nrow(priors)), function(i) {
  data.frame(
    p = p_grid,
    density = dbeta(p_grid, priors$a[i] + y_small, priors$b[i] + n_small - y_small),
    prior = priors$prior[i]
  )
}))

ggplot(plot_prior_sens, aes(p, density, color = prior)) +
  geom_line(linewidth = 1.1) +
  labs(
    title = "Prior sensitivity when data are small",
    subtitle = "Different priors lead to visibly different posteriors when n is small.",
    x = "p", y = "posterior density", color = "prior"
  )

12 Bayesian decision and loss

Bayesian inference becomes decision-making when we specify an action \(a\) and a loss function \(L(a,\theta)\). The Bayes action minimizes posterior expected loss:

\[ a^*(y)=\arg\min_a E\{L(a,\theta)\mid y\}. \]

Examples:

  • squared error loss \(L(a,\theta)=(a-\theta)^2\) gives posterior mean;
  • absolute error loss \(L(a,\theta)=|a-\theta|\) gives posterior median;
  • asymmetric loss can be used when underestimating a hazard is worse than overestimating it.

In natural science, decision problems include observing-time allocation, alert thresholds, telescope follow-up scheduling, quality control, hazard ranking, and experiment design.

13 Compact Bayesian workflow for research

A responsible Bayesian analysis usually follows this loop:

  1. Scientific question: What is the estimand or latent state?
  2. Probability model: What is random, what is observed, what is latent?
  3. Prior specification: What information or regularization is reasonable?
  4. Posterior computation: Conjugacy, grid, MCMC, variational Bayes, or sequential Monte Carlo.
  5. Posterior summaries: estimates, credible intervals, posterior probabilities.
  6. Posterior predictive checking: can the fitted model reproduce important features of the data?
  7. Sensitivity analysis: do conclusions change under defensible prior/model changes?
  8. Scientific interpretation: translate probability statements back to the domain problem.

Bayesian inference is not the claim that a model is true. It is a coherent way to reason conditionally on a model. Scientific honesty requires checking and revising the model.

14 Session information

sessionInfo()
## R version 4.4.0 (2024-04-24 ucrt)
## Platform: x86_64-w64-mingw32/x64
## Running under: Windows 11 x64 (build 26200)
## 
## Matrix products: default
## 
## 
## locale:
## [1] LC_COLLATE=English_United States.utf8 
## [2] LC_CTYPE=English_United States.utf8   
## [3] LC_MONETARY=English_United States.utf8
## [4] LC_NUMERIC=C                          
## [5] LC_TIME=English_United States.utf8    
## 
## time zone: Asia/Calcutta
## tzcode source: internal
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
## [1] MASS_7.3-65   knitr_1.51    ggplot2_4.0.2
## 
## loaded via a namespace (and not attached):
##  [1] vctrs_0.6.5        cli_3.6.2          rlang_1.1.4        xfun_0.56         
##  [5] otel_0.2.0         generics_0.1.4     S7_0.2.0           jsonlite_1.8.8    
##  [9] labeling_0.4.3     glue_1.7.0         htmltools_0.5.9    sass_0.4.10       
## [13] scales_1.4.0       rmarkdown_2.30     grid_4.4.0         tibble_3.2.1      
## [17] evaluate_1.0.5     jquerylib_0.1.4    fastmap_1.2.0      yaml_2.3.12       
## [21] lifecycle_1.0.5    compiler_4.4.0     dplyr_1.1.4        RColorBrewer_1.1-3
## [25] pkgconfig_2.0.3    rstudioapi_0.18.0  farver_2.1.2       digest_0.6.35     
## [29] R6_2.6.1           tidyselect_1.2.1   dichromat_2.0-0.1  pillar_1.11.1     
## [33] magrittr_2.0.3     bslib_0.10.0       withr_3.0.2        tools_4.4.0       
## [37] gtable_0.3.6       cachem_1.1.0

15 References and further reading

Berger, J. O. (1985). Statistical Decision Theory and Bayesian Analysis. Springer.

Gelman, A., Carlin, J. B., Stern, H. S., Dunson, D. B., Vehtari, A., and Rubin, D. B. (2013). Bayesian Data Analysis. Chapman and Hall/CRC.

Gelman, A., Vehtari, A., Simpson, D., et al. (2020). Bayesian workflow. arXiv:2011.01808.

Gelman, A. and Shalizi, C. R. (2013). Philosophy and the practice of Bayesian statistics. British Journal of Mathematical and Statistical Psychology, 66, 8–38.

Jaynes, E. T. (2003). Probability Theory: The Logic of Science. Cambridge University Press.

Robert, C. P. (2007). The Bayesian Choice. Springer.

Venables, W. N. and Ripley, B. D. (2002). Modern Applied Statistics with S. Springer.