required_pkgs <- c("ggplot2", "dplyr", "tidyr", "tibble", "knitr", "purrr", "scales")
to_install <- setdiff(required_pkgs, rownames(installed.packages()))
if(length(to_install) > 0) install.packages(to_install, repos = "https://cloud.r-project.org")

suppressPackageStartupMessages({
  library(ggplot2)
  library(dplyr)
  library(tidyr)
  library(tibble)
  library(knitr)
  library(purrr)
  library(scales)
})

set.seed(20260708)
theme_set(theme_minimal(base_size = 13))

1 Why this note?

This note is a foundation-first continuation tutorial on statistical inference and data science for basic sciences.

The goal is to explain why probability is not merely a chapter before statistics, but the mathematical language behind:

The intended audience is a science researcher or M.Sc./Ph.D. student who already knows elementary probability, but wants to understand the backbone more deeply.

Main philosophy:
In natural science, data are not just tables. Data are partial observations of uncertain processes. Probability tells us what can be measured, distributions tell us what uncertainty looks like, entropy tells us how much uncertainty remains, and Bayesian inference tells us how uncertainty changes after observing data.

2 Measure-theoretic probability: the backbone

2.1 Probability space

A modern probability model begins with a triple

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

Here:

  • \(\Omega\) is the sample space: the set of all possible elementary outcomes.
  • \(\mathcal F\) is a sigma-algebra: the collection of events to which probabilities can be assigned.
  • \(P\) is a probability measure: a function assigning probabilities to events.

The measure \(P\) satisfies:

\[ P(\Omega)=1, \]

and if \(A_1,A_2,\ldots\) are disjoint events,

\[ P\left(\bigcup_{i=1}^{\infty}A_i\right)=\sum_{i=1}^{\infty}P(A_i). \]

This countable additivity is not a technical decoration. It is what makes limits, stochastic processes, convergence, and continuous data mathematically stable.

2.2 Why do we need sigma-algebras?

In finite probability, we can assign probability to every subset of \(\Omega\). But in continuous sample spaces, assigning probabilities to every possible subset leads to paradoxes and non-measurable sets.

Therefore, we specify a scientifically meaningful class of events \(\mathcal F\). For real-valued observations, the usual choice is the Borel sigma-algebra, generated by intervals.

Scientific translation:
We only assign probabilities to events that are observable or mathematically well-defined.

2.3 Random variable as a measurable function

A random variable is not originally a variable in the ordinary algebraic sense. It is a measurable map

\[ X:(\Omega,\mathcal F)\rightarrow (\mathbb R,\mathcal B(\mathbb R)). \]

For every Borel set \(B\subseteq \mathbb R\), the event

\[ \{\omega:X(\omega)\in B\} \]

must belong to \(\mathcal F\).

Example:
If \(X\) is tomorrow’s temperature, then \(\{X>35\}\) must be an event whose probability is meaningful.

2.4 Distribution as push-forward measure

The distribution of \(X\) is the probability measure \(P_X\) defined by

\[ P_X(B)=P(X\in B)=P\{\omega:X(\omega)\in B\}. \]

This is called the push-forward of \(P\) under \(X\).

So a distribution is not just a curve. It is a measure induced by a random variable.

2.5 Expectation as an integral

The expectation is an integral with respect to a probability measure:

\[ \mathbb E[X]=\int_{\Omega}X(\omega)dP(\omega). \]

Equivalently, if \(X\) has distribution \(P_X\),

\[ \mathbb E[g(X)]=\int_{\mathbb R}g(x)dP_X(x). \]

This formula is the basis of almost everything in inference: moments, loss, risk, likelihood, posterior prediction, filtering, and decision theory.

3 Probability language for scientific data

A scientific measurement can be written as

\[ Y = \text{signal} + \text{noise}. \]

But in probability language, we say:

\[ Y\sim P_\theta, \]

where \(P_\theta\) is a distribution indexed by unknown scientific parameter \(\theta\).

Examples:

Domain Observation Possible random variable
Physics detector count Poisson or negative binomial count
Astronomy photon arrival time Poisson process waiting time
Geoscience earthquake magnitude heavy-tailed or threshold distribution
Chemistry reaction waiting time exponential/gamma-like variable
Biology cell count or expression count distribution with overdispersion
Materials measured strength normal/lognormal/extreme-value variable

A distribution should not be chosen only because it is familiar. It should reflect:

4 Entropy: what exactly does it mean?

Entropy is a quantitative measure of uncertainty or spread of a probability distribution.

4.1 Discrete Shannon entropy

For a discrete distribution \(p=(p_1,\ldots,p_K)\),

\[ H(p)=-\sum_{k=1}^{K}p_k\log p_k. \]

If the logarithm is natural, entropy is measured in nats. If the logarithm is base 2, entropy is measured in bits.

Interpretation:

  • Low entropy: distribution is concentrated; uncertainty is small.
  • High entropy: distribution is spread out; uncertainty is large.
  • Maximum entropy on a finite set occurs at the uniform distribution.
entropy_disc <- function(p, base = exp(1)) {
  p <- p[p > 0]
  -sum(p * log(p)) / log(base)
}

diff_entropy_grid <- function(x, f) {
  dx <- mean(diff(x))
  f <- pmax(f, 0)
  f <- f / sum(f * dx)
  positive <- f > 0
  -sum(f[positive] * log(f[positive])) * dx
}

4.2 Example: entropy of a Bernoulli event

Suppose a scientific event has probability \(p\). For instance, a detector flags a rare event, an earthquake exceeds a threshold, or a patient responds to treatment.

The entropy is

\[ H(p)=-p\log_2p-(1-p)\log_2(1-p). \]

p_grid <- seq(0.001, 0.999, length.out = 500)
bern_entropy <- tibble::tibble(
  p = p_grid,
  entropy_bits = -p * log2(p) - (1 - p) * log2(1 - p)
)

ggplot(bern_entropy, aes(p, entropy_bits)) +
  geom_line(linewidth = 1.1) +
  geom_vline(xintercept = 0.5, linetype = "dashed") +
  labs(
    title = "Entropy of a Bernoulli event",
    subtitle = "Maximum uncertainty occurs at p = 0.5, not at rare or almost certain events",
    x = "Probability of event",
    y = "Entropy in bits"
  )

Scientific message:
An event with probability 0.5 is maximally uncertain. A rare event may be scientifically important, but its one-step Shannon entropy is not maximal.

4.3 Continuous differential entropy

For a continuous density \(f\), differential entropy is

\[ h(f)=-\int f(x)\log f(x)dx. \]

Important caution:

  • Discrete entropy is invariant to relabelling.
  • Differential entropy depends on the coordinate scale.
  • Therefore, for continuous problems, a more principled formulation uses relative entropy with respect to a reference measure or baseline density.

4.4 Relative entropy / Kullback–Leibler divergence

For two densities \(p\) and \(m\),

\[ D_{\mathrm{KL}}(p\|m)=\int p(x)\log\frac{p(x)}{m(x)}dx. \]

This measures how far \(p\) is from reference density \(m\). In maximum-entropy reasoning, the continuous invariant form often maximizes

\[ -\int p(x)\log\frac{p(x)}{m(x)}dx \]

subject to known constraints.

Scientific message:
Entropy is not magic randomness. It is a rule for choosing the least committed distribution consistent with known information.

5 Maximum entropy: what distribution follows from what information?

The maximum-entropy principle says:

Among all distributions satisfying the information we truly know, choose the one with largest entropy.

This does not mean nature always maximizes entropy in every statistical model. It means that if our information is limited to certain constraints, the maximum-entropy distribution is the least extra-assumptive distribution compatible with those constraints.

5.1 Core maximum-entropy table

maxent_table <- tibble::tribble(
  ~Support, ~Known_constraints, ~Maximum_entropy_distribution, ~Natural_science_interpretation,
  "Finite set {1,...,K}", "Only probabilities sum to 1", "Uniform", "No state is preferred before observing more information",
  "[a,b]", "Only bounded support", "Uniform(a,b)", "Measurement known only to lie inside physical range",
  "[0,infinity)", "Fixed mean E[X] = mu", "Exponential(rate = 1/mu)", "Waiting time with known average and no further structure",
  "R", "Fixed mean and variance", "Normal(mu, sigma^2)", "Measurement error with known first two moments",
  "R^d", "Fixed mean vector and covariance matrix", "Multivariate normal", "Correlated measurement errors or multivariate fluctuations",
  "Nonnegative integers", "Fixed mean", "Geometric", "Discrete waiting/count until first success under mean constraint",
  "Finite energy states", "Fixed expected energy", "Gibbs/Boltzmann", "Equilibrium statistical mechanics style uncertainty over states",
  "Positive real line", "Fixed E[X] and E[log X]", "Gamma", "Positive measurements with scale and shape information",
  "[0,1]", "Fixed E[log X] and E[log(1-X)]", "Beta", "Uncertain proportion with boundary-sensitive information",
  "Circle", "Fixed mean resultant direction", "von Mises", "Angular data with preferred direction and concentration"
)
knitr::kable(maxent_table, caption = "Maximum-entropy distributions and their information constraints")
Maximum-entropy distributions and their information constraints
Support Known_constraints Maximum_entropy_distribution Natural_science_interpretation
Finite set {1,…,K} Only probabilities sum to 1 Uniform No state is preferred before observing more information
[a,b] Only bounded support Uniform(a,b) Measurement known only to lie inside physical range
[0,infinity) Fixed mean E[X] = mu Exponential(rate = 1/mu) Waiting time with known average and no further structure
R Fixed mean and variance Normal(mu, sigma^2) Measurement error with known first two moments
R^d Fixed mean vector and covariance matrix Multivariate normal Correlated measurement errors or multivariate fluctuations
Nonnegative integers Fixed mean Geometric Discrete waiting/count until first success under mean constraint
Finite energy states Fixed expected energy Gibbs/Boltzmann Equilibrium statistical mechanics style uncertainty over states
Positive real line Fixed E[X] and E[log X] Gamma Positive measurements with scale and shape information
[0,1] Fixed E[log X] and E[log(1-X)] Beta Uncertain proportion with boundary-sensitive information
Circle Fixed mean resultant direction von Mises Angular data with preferred direction and concentration

6 Finite maximum entropy: uniform distribution

If the only information is that one of \(K\) states occurs, then the maximum entropy distribution is uniform:

\[ p_1=\cdots=p_K=1/K. \]

K <- 6
random_probs <- replicate(800, {
  z <- rexp(K)
  z / sum(z)
})
H_random <- apply(random_probs, 2, entropy_disc)
H_uniform <- entropy_disc(rep(1 / K, K))

entropy_df <- tibble::tibble(entropy = H_random)

ggplot(entropy_df, aes(entropy)) +
  geom_histogram(bins = 35, fill = "grey80", color = "white") +
  geom_vline(xintercept = H_uniform, linewidth = 1.2, linetype = "dashed") +
  labs(
    title = "Random discrete distributions have entropy below the uniform maximum",
    subtitle = paste0("For K = ", K, ", uniform entropy = log(K) = ", round(H_uniform, 3), " nats"),
    x = "Entropy in nats",
    y = "Number of random distributions"
  )

7 Normal distribution as maximum entropy under mean and variance

If we only know

\[ \mathbb E[X]=\mu,\qquad \mathrm{Var}(X)=\sigma^2, \]

then the maximum differential entropy distribution on \(\mathbb R\) is

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

This is why normal errors are not merely convenient. They are the least informative continuous error model when only mean and variance are specified.

7.1 Numerical illustration: same mean and variance, different entropy

x <- seq(-6, 6, length.out = 5000)

# All distributions below are centered around 0 and scaled approximately to variance 1.
normal_d <- dnorm(x, mean = 0, sd = 1)
uniform_d <- dunif(x, min = -sqrt(3), max = sqrt(3))

laplace_density <- function(x, b) exp(-abs(x) / b) / (2 * b)
laplace_d <- laplace_density(x, b = 1 / sqrt(2))  # variance = 1

# Student t with df=5 scaled to variance 1.
df_t <- 5
scale_t <- sqrt((df_t - 2) / df_t)
t_d <- dt(x / scale_t, df = df_t) / scale_t

density_df <- tibble::tibble(
  x = x,
  Normal = normal_d,
  Uniform = uniform_d,
  Laplace = laplace_d,
  Student_t5 = t_d
) |>
  tidyr::pivot_longer(-x, names_to = "distribution", values_to = "density")

entropy_compare <- tibble::tibble(
  distribution = c("Normal", "Uniform", "Laplace", "Student_t5_scaled"),
  entropy_nats = c(
    diff_entropy_grid(x, normal_d),
    diff_entropy_grid(x, uniform_d),
    diff_entropy_grid(x, laplace_d),
    diff_entropy_grid(x, t_d)
  )
) |>
  dplyr::arrange(dplyr::desc(entropy_nats))

knitr::kable(entropy_compare, digits = 4, caption = "Approximate differential entropy among distributions scaled to comparable variance")
Approximate differential entropy among distributions scaled to comparable variance
distribution entropy_nats
Normal 1.4189
Student_t5_scaled 1.3668
Laplace 1.3446
Uniform 1.2431 
ggplot(density_df, aes(x, density, color = distribution)) +
  geom_line(linewidth = 1) +
  labs(
    title = "Same broad scale, different tail assumptions",
    subtitle = "The normal is the maximum-entropy model when mean and variance are the only constraints on the real line",
    x = "x",
    y = "Density"
  )

Scientific message:
If you assume normal errors, you are implicitly saying: beyond mean and variance, I do not encode further information. But if the science says heavy tails, bounded support, or asymmetry, a different distribution may be more appropriate.

8 Exponential distribution as maximum entropy on positive support

If \(X\ge 0\) and only the mean \(\mathbb E[X]=\mu\) is known, then the maximum-entropy distribution is exponential:

\[ f(x)=\frac{1}{\mu}\exp(-x/\mu),\qquad x\ge 0. \]

This is important for waiting times: radioactive decay waiting times, reaction waiting times, arrival intervals, and simple lifetime models.

x_pos <- seq(0, 8, length.out = 4000)
mu <- 1

exp_d <- dexp(x_pos, rate = 1 / mu)
# Gamma with mean 1, shape 3, scale 1/3
gamma_d <- dgamma(x_pos, shape = 3, scale = 1 / 3)
# Half-normal with mean 1
sigma_hn <- sqrt(pi / 2)
halfnorm_d <- sqrt(2) / (sigma_hn * sqrt(pi)) * exp(-x_pos^2 / (2 * sigma_hn^2))
# Uniform on [0,2] has mean 1
unif_pos_d <- dunif(x_pos, 0, 2)

positive_entropy <- tibble::tibble(
  distribution = c("Exponential", "Gamma_shape3", "Half_normal", "Uniform_0_2"),
  approx_entropy_nats = c(
    diff_entropy_grid(x_pos, exp_d),
    diff_entropy_grid(x_pos, gamma_d),
    diff_entropy_grid(x_pos, halfnorm_d),
    diff_entropy_grid(x_pos, unif_pos_d)
  )
) |>
  dplyr::arrange(dplyr::desc(approx_entropy_nats))

knitr::kable(positive_entropy, digits = 4, caption = "Positive distributions with mean about 1: exponential has maximum entropy under only the mean constraint")
Positive distributions with mean about 1: exponential has maximum entropy under only the mean constraint
distribution approx_entropy_nats
Exponential 0.9970
Half_normal 0.9519
Gamma_shape3 0.7490
Uniform_0_2 0.6934 
positive_df <- tibble::tibble(
  x = x_pos,
  Exponential = exp_d,
  Gamma_shape3 = gamma_d,
  Half_normal = halfnorm_d,
  Uniform_0_2 = unif_pos_d
) |>
  tidyr::pivot_longer(-x, names_to = "distribution", values_to = "density")

ggplot(positive_df, aes(x, density, color = distribution)) +
  geom_line(linewidth = 1) +
  labs(
    title = "Positive waiting-time models with the same mean tell different scientific stories",
    subtitle = "Exponential is the maximum-entropy choice if only positivity and mean are known",
    x = "Waiting time",
    y = "Density"
  )

9 Gibbs/Boltzmann distribution from maximum entropy

Suppose a system can occupy finite states \(i=1,\ldots,K\) with energies \(E_i\). If we know only the expected energy,

\[ \sum_i p_iE_i=\bar E, \]

then maximizing entropy gives

\[ p_i=\frac{\exp(-\beta E_i)}{Z(\beta)}, \qquad Z(\beta)=\sum_j\exp(-\beta E_j). \]

This is the Gibbs/Boltzmann form. The parameter \(\beta\) is chosen to match the expected-energy constraint.

energy <- 0:7
betas <- c(0, 0.3, 0.8, 1.5)

boltzmann <- function(beta, E) {
  p <- exp(-beta * E)
  p / sum(p)
}

boltz_df <- purrr::map_dfr(betas, function(b) {
  p <- boltzmann(b, energy)
  tibble::tibble(beta = b, energy = energy, probability = p)
})

boltz_summary <- boltz_df |>
  dplyr::group_by(beta) |>
  dplyr::summarise(
    mean_energy = sum(energy * probability),
    entropy_nats = entropy_disc(probability),
    .groups = "drop"
  )

knitr::kable(boltz_summary, digits = 4, caption = "Boltzmann distributions: increasing beta concentrates probability on low-energy states")
Boltzmann distributions: increasing beta concentrates probability on low-energy states
beta mean_energy entropy_nats
0.0 3.5000 2.0794
0.3 2.0601 1.8732
0.8 0.8027 1.2371
1.5 0.2872 0.6832 
ggplot(boltz_df, aes(factor(energy), probability, fill = factor(beta))) +
  geom_col(position = "dodge") +
  labs(
    title = "Maximum entropy with expected-energy constraint",
    subtitle = "Higher beta means stronger concentration on lower-energy states",
    x = "Energy level",
    y = "Probability",
    fill = "beta"
  )

Scientific message:
The same entropy idea connects statistics with statistical mechanics. A distribution is determined not by arbitrary curve fitting, but by the information constraints we choose to honour.

10 Circular maximum entropy: von Mises distribution

Natural science often has angular observations:

On the circle, a natural maximum-entropy distribution under fixed mean resultant direction is the von Mises distribution:

\[ f(\theta\mid \mu,\kappa)=\frac{1}{2\pi I_0(\kappa)}\exp\{\kappa\cos(\theta-\mu)\}. \]

Here:

dvonmises <- function(theta, mu = 0, kappa = 1) {
  exp(kappa * cos(theta - mu)) / (2 * pi * besselI(kappa, nu = 0))
}

theta <- seq(-pi, pi, length.out = 1000)
vm_df <- tibble::tibble(
  theta = theta,
  kappa_0 = dvonmises(theta, 0, 0),
  kappa_1 = dvonmises(theta, 0, 1),
  kappa_4 = dvonmises(theta, 0, 4)
) |>
  tidyr::pivot_longer(-theta, names_to = "model", values_to = "density")

ggplot(vm_df, aes(theta, density, color = model)) +
  geom_line(linewidth = 1) +
  scale_x_continuous(labels = function(z) round(z / pi, 2), breaks = seq(-pi, pi, by = pi/2)) +
  labs(
    title = "von Mises family on the circle",
    subtitle = "kappa = 0 is uniform; larger kappa means stronger directional concentration",
    x = "Angle in multiples of pi",
    y = "Density"
  )

11 Real natural-science example 1: precipitation and distribution choice

The built-in precip data contain average annual precipitation values for many U.S. cities. It is a simple positive-valued environmental dataset.

The purpose here is not to build a final climate model. It is to show how distribution choice reflects scientific assumptions:

data("precip")
prec <- as.numeric(precip)
prec_df <- tibble::tibble(precip_inch = prec)

mean_prec <- mean(prec)
var_prec <- var(prec)

# Moment-calibrated gamma.
gamma_shape <- mean_prec^2 / var_prec
gamma_scale <- var_prec / mean_prec

# Densities for plotting.
xp <- seq(min(prec) * 0.7, max(prec) * 1.2, length.out = 800)
model_dens <- tibble::tibble(
  x = xp,
  Normal = dnorm(xp, mean = mean_prec, sd = sd(prec)),
  Gamma = dgamma(xp, shape = gamma_shape, scale = gamma_scale),
  Exponential = dexp(xp, rate = 1 / mean_prec)
) |>
  tidyr::pivot_longer(-x, names_to = "model", values_to = "density")

# Log-likelihood and AIC-like comparison.
ll_normal <- sum(dnorm(prec, mean_prec, sd(prec), log = TRUE))
ll_gamma <- sum(dgamma(prec, shape = gamma_shape, scale = gamma_scale, log = TRUE))
ll_exp <- sum(dexp(prec, rate = 1 / mean_prec, log = TRUE))

model_table <- tibble::tibble(
  model = c("Normal(moment)", "Gamma(moment)", "Exponential(mean)"),
  parameters_used = c(2, 2, 1),
  log_likelihood = c(ll_normal, ll_gamma, ll_exp),
  AIC_like = -2 * log_likelihood + 2 * parameters_used
) |>
  dplyr::arrange(AIC_like)

knitr::kable(model_table, digits = 3, caption = "Simple distribution comparison for positive precipitation data")
Simple distribution comparison for positive precipitation data
model parameters_used log_likelihood AIC_like
Normal(moment) 2 -282.077 568.155
Gamma(moment) 2 -290.552 585.104
Exponential(mean) 1 -318.645 639.291 
ggplot(prec_df, aes(precip_inch)) +
  geom_histogram(aes(y = after_stat(density)), bins = 18, fill = "grey85", color = "white") +
  geom_line(data = model_dens, aes(x, density, color = model), linewidth = 1.1) +
  labs(
    title = "Distribution choice in environmental data: precipitation is positive and skewed",
    subtitle = "Normal is not automatically appropriate; support and skewness matter",
    x = "Average annual precipitation",
    y = "Density"
  )

Motivation:
Maximum entropy helps us ask: what information do we truly know? Positivity alone plus mean suggests exponential. Positivity plus additional shape information may suggest gamma. Symmetric errors with mean and variance suggest normal.

12 Stochastic process: random variables indexed by time

A stochastic process is a collection

\[ \{X_t:t\in T\}, \]

where each \(X_t\) is a random variable.

Examples:

12.1 Filtration: information grows over time

A filtration is an increasing family of sigma-algebras

\[ \mathcal F_0\subseteq \mathcal F_1\subseteq\cdots\subseteq \mathcal F_t\subseteq\cdots. \]

\(\mathcal F_t\) represents all information available up to time \(t\).

A process \(X_t\) is adapted to \(\mathcal F_t\) if \(X_t\) is known at time \(t\).

Scientific meaning:
A filtering algorithm should not use future observations to estimate the present state. It must respect the filtration.

13 Real natural-science example 2: ecological stochastic process with lynx data

The built-in lynx data contain annual Canadian lynx trappings. It is a classic ecological time-series dataset with cyclic behaviour.

data("lynx")
lynx_df <- tibble::tibble(
  year = as.numeric(time(lynx)),
  count = as.numeric(lynx),
  log_count = log(as.numeric(lynx))
)

ggplot(lynx_df, aes(year, count)) +
  geom_line(linewidth = 0.9) +
  geom_point(size = 1.2) +
  scale_y_log10(labels = scales::comma) +
  labs(
    title = "Canadian lynx trappings as an ecological stochastic process",
    subtitle = "The process is not independent noise; it has temporal structure and cycles",
    x = "Year",
    y = "Count on log scale"
  )

lynx_ar <- stats::arima(log(lynx), order = c(2, 0, 0))
lynx_fitted <- as.numeric(log(lynx) - residuals(lynx_ar))
lynx_fit_df <- lynx_df |>
  dplyr::mutate(fitted_log_count = lynx_fitted)

ar_table <- tibble::tibble(
  quantity = c("AR coefficient 1", "AR coefficient 2", "intercept/mean", "innovation variance"),
  estimate = c(lynx_ar$coef[1], lynx_ar$coef[2], lynx_ar$coef[3], lynx_ar$sigma2)
)
knitr::kable(ar_table, digits = 4, caption = "Simple AR(2) model fitted to log lynx counts")
Simple AR(2) model fitted to log lynx counts
quantity estimate
AR coefficient 1 1.3776
AR coefficient 2 -0.7399
intercept/mean 6.6863
innovation variance 0.2708 
ggplot(lynx_fit_df, aes(year)) +
  geom_line(aes(y = log_count, color = "Observed log count"), linewidth = 0.9) +
  geom_line(aes(y = fitted_log_count, color = "AR(2) fitted value"), linewidth = 0.9) +
  labs(
    title = "A basic stochastic-process model: AR(2) structure",
    subtitle = "Even a simple time-series model captures part of the cyclic dependence",
    x = "Year",
    y = "log(count)",
    color = "Series"
  )

Scientific message:
When data are temporally dependent, treating observations as independent destroys part of the scientific structure. Stochastic process models encode dependence.

14 Filtering as repeated probabilistic inference

Filtering means estimating the current hidden state using observations available up to the current time:

\[ p(x_t\mid y_1,\ldots,y_t). \]

This is a Bayesian object even before we call it Bayesian statistics. Filtering is probabilistic conditioning over time.

14.1 A nonlinear particle-filter example

Suppose a hidden scientific state evolves as

\[ X_t=0.7X_{t-1}+\epsilon_t, \]

but the sensor observes

\[ Y_t=\sin(X_t)+\eta_t. \]

This is nonlinear. A Kalman filter is not directly exact. A particle filter approximates the filtering distribution by many simulated particles.

wtd_quantile <- function(x, w, probs = c(0.05, 0.5, 0.95)) {
  ord <- order(x)
  x <- x[ord]
  w <- w[ord] / sum(w)
  cw <- cumsum(w)
  approx(cw, x, xout = probs, rule = 2, ties = "ordered")$y
}
set.seed(77)
Tn <- 90
sigma_state <- 0.45
sigma_obs <- 0.25

x_true <- numeric(Tn)
y_obs <- numeric(Tn)
x_true[1] <- rnorm(1, 0, 1)
y_obs[1] <- sin(x_true[1]) + rnorm(1, 0, sigma_obs)

for(t in 2:Tn) {
  x_true[t] <- 0.7 * x_true[t - 1] + rnorm(1, 0, sigma_state)
  y_obs[t] <- sin(x_true[t]) + rnorm(1, 0, sigma_obs)
}

N <- 1500
particles <- rnorm(N, 0, 1)
weights <- rep(1 / N, N)

pf_store <- vector("list", Tn)

for(t in 1:Tn) {
  if(t > 1) particles <- 0.7 * particles + rnorm(N, 0, sigma_state)
  
  logw <- dnorm(y_obs[t], mean = sin(particles), sd = sigma_obs, log = TRUE)
  logw <- logw - max(logw)
  weights <- exp(logw)
  weights <- weights / sum(weights)
  ess <- 1 / sum(weights^2)
  qs <- wtd_quantile(particles, weights, probs = c(0.05, 0.5, 0.95))
  
  pf_store[[t]] <- tibble::tibble(
    time = t,
    true_state = x_true[t],
    observation = y_obs[t],
    estimate = sum(weights * particles),
    q05 = qs[1],
    q50 = qs[2],
    q95 = qs[3],
    ess = ess
  )
  
  if(ess < N / 2) {
    idx <- sample.int(N, size = N, replace = TRUE, prob = weights)
    particles <- particles[idx]
    weights <- rep(1 / N, N)
  }
}

pf_df <- dplyr::bind_rows(pf_store)

knitr::kable(head(pf_df, 8), digits = 3, caption = "First few filtering summaries from the particle filter")
First few filtering summaries from the particle filter
time true_state observation estimate q05 q50 q95 ess
1 -0.550 -0.250 -0.276 -0.779 -0.253 0.159 562.046
2 -0.097 0.164 0.087 -0.278 0.074 0.477 780.039
3 0.009 0.293 0.242 -0.147 0.232 0.671 685.838
4 -0.431 -0.451 -0.325 -0.707 -0.337 0.046 516.395
5 -0.236 0.127 0.044 -0.324 0.047 0.423 784.129
6 -1.489 -1.057 -1.014 -1.593 -0.994 -0.519 283.548
7 -1.105 -0.902 -1.026 -1.628 -0.993 -0.534 986.034
8 -0.648 -0.456 -0.541 -1.061 -0.529 -0.087 960.435 
ggplot(pf_df, aes(time)) +
  geom_ribbon(aes(ymin = q05, ymax = q95), fill = "grey80", alpha = 0.7) +
  geom_line(aes(y = true_state, color = "True hidden state"), linewidth = 1) +
  geom_line(aes(y = estimate, color = "Particle-filter estimate"), linewidth = 1) +
  geom_point(aes(y = observation, color = "Noisy observation"), size = 1.2, alpha = 0.75) +
  labs(
    title = "Filtering is sequential probabilistic inference",
    subtitle = "The grey band is a filtering uncertainty band for the hidden state",
    x = "Time",
    y = "State / observation scale",
    color = "Quantity"
  )

Scientific message:
Filtering is not only an engineering algorithm. It is a probability statement evolving with time:

\[ p(x_t\mid \mathcal F_t^Y), \]

where \(\mathcal F_t^Y\) is the information generated by observations \(Y_1,\ldots,Y_t\).

15 Maximum entropy and model choice: a warning

Maximum entropy does not say that the chosen distribution is always true. It says that among all distributions satisfying the specified constraints, it is least informative.

If the constraints are wrong, the model can be wrong.

Examples:

Therefore, scientific modelling is an iterative workflow:

  1. Identify support and constraints.
  2. Choose a probabilistic model.
  3. Fit parameters.
  4. Check residuals and predictions.
  5. Revise constraints or model family.

16 Philosophy of Bayesian inference

Bayesian inference is often called probabilistic inference because it uses probability to represent uncertainty about unknown quantities.

Bayes’ theorem is

\[ p(\theta\mid y)=\frac{p(y\mid \theta)p(\theta)}{p(y)}. \]

Here:

16.1 Subjective Bayes philosophy

Subjective Bayes interprets probability as a rational degree of belief.

This does not mean arbitrary opinion. A well-formed subjective prior should be:

  • coherent,
  • transparent,
  • based on available scientific knowledge,
  • and checked against data.

Example:
If an astronomer has prior information about expected background rate, or a geoscientist has expert knowledge about plausible fault-slip parameters, that information can enter the prior.

Subjective Bayes is powerful when expert knowledge is real and valuable.

16.2 Objective Bayes philosophy

Objective Bayes attempts to choose priors by formal rules, such as:

  • symmetry,
  • invariance,
  • reference priors,
  • Jeffreys priors,
  • maximum-entropy priors under stated constraints.

Example:
For a positive scale parameter \(\sigma\), a common objective-style prior is proportional to \(1/\sigma\), reflecting scale invariance.

Objective Bayes is useful when subjective prior information is weak or when a default analysis is desired.

16.3 Practical Bayes: modelling plus checking

In modern scientific data analysis, Bayesian inference is not only prior times likelihood. It is a workflow:

\[ \text{prior}\rightarrow\text{likelihood}\rightarrow\text{posterior}\rightarrow\text{posterior predictive check}\rightarrow\text{model revision}. \]

This is important because all models are approximations.

17 Tiny Bayesian demonstration: coin probability with different priors

Suppose we observe \(y\) successes in \(n\) trials. With a beta prior

\[ \theta\sim \mathrm{Beta}(a,b), \]

and binomial likelihood

\[ Y\mid\theta\sim\mathrm{Binomial}(n,\theta), \]

posterior is

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

n <- 20
y <- 14
prior_grid <- tibble::tribble(
  ~prior, ~a, ~b,
  "Uniform Beta(1,1)", 1, 1,
  "Jeffreys Beta(1/2,1/2)", 0.5, 0.5,
  "Skeptical Beta(8,8)", 8, 8
)

theta_grid <- seq(0.001, 0.999, length.out = 1000)

posterior_df <- prior_grid |>
  dplyr::mutate(post_a = a + y, post_b = b + n - y) |>
  tidyr::crossing(theta = theta_grid) |>
  dplyr::mutate(density = dbeta(theta, post_a, post_b))

posterior_summary <- prior_grid |>
  dplyr::mutate(
    post_a = a + y,
    post_b = b + n - y,
    posterior_mean = post_a / (post_a + post_b),
    q025 = qbeta(0.025, post_a, post_b),
    q975 = qbeta(0.975, post_a, post_b),
    prob_theta_gt_half = 1 - pbeta(0.5, post_a, post_b)
  ) |>
  dplyr::select(prior, posterior_mean, q025, q975, prob_theta_gt_half)

knitr::kable(posterior_summary, digits = 3, caption = "Posterior summaries after observing 14 successes in 20 trials")
Posterior summaries after observing 14 successes in 20 trials
prior posterior_mean q025 q975 prob_theta_gt_half
Uniform Beta(1,1) 0.682 0.478 0.854 0.961
Jeffreys Beta(1/2,1/2) 0.690 0.483 0.864 0.965
Skeptical Beta(8,8) 0.611 0.449 0.761 0.912 
ggplot(posterior_df, aes(theta, density, color = prior)) +
  geom_line(linewidth = 1.1) +
  geom_vline(xintercept = y / n, linetype = "dashed") +
  labs(
    title = "Bayesian updating: same data, different transparent priors",
    subtitle = "With more data, reasonable priors typically become less influential",
    x = "Unknown probability theta",
    y = "Posterior density"
  )

Interpretation:
Bayesian inference does not hide assumptions. It makes them explicit. Different priors can be compared, justified, and stress-tested.

18 Putting the philosophy together

18.1 Frequentist, ML, maximum entropy, and Bayes

Framework Primary question Strength Risk if misused
Frequentist inference What would happen under repeated sampling? Calibration, tests, confidence procedures Can become mechanical if model assumptions ignored
Machine learning What predicts well on future data? Flexible prediction and pattern recognition Can exploit shortcuts and fail under shift
Maximum entropy What is the least committed distribution under constraints? Principled distribution construction Wrong constraints give wrong model
Bayesian inference How should uncertainty update after data? Full uncertainty propagation and decision support Priors and likelihood must be checked

A mature scientific workflow uses all of them wisely.

19 Suggested discussion questions for students

  1. If data are positive and skewed, why might a normal model be scientifically questionable?
  2. What constraint leads to the normal distribution by maximum entropy?
  3. Why is differential entropy more delicate than discrete entropy?
  4. What does a filtration represent in sequential experiments?
  5. When should expert knowledge enter a prior?
  6. What is the difference between subjective prior and arbitrary prior?
  7. Why does objective Bayes often rely on invariance?
  8. How would you check whether a maximum-entropy model is adequate?

20 Final message

Measure-theoretic probability gives the mathematical foundation. Entropy explains uncertainty and principled distribution construction. Maximum entropy tells us how common distributions arise from constraints. Stochastic processes and filtrations explain evolving information. Bayesian inference then becomes a natural next step: uncertainty is updated as data arrive.

For basic sciences, the message is simple but powerful:

Do not merely fit a curve. State the sample space, identify the measurable quantity, choose distributions from scientific constraints, quantify uncertainty, and update that uncertainty when new evidence arrives.