Probability and Statistics: A Geometric Foundation

1. Probability as measure

Most first encounters with probability present it as a set of rules for computing with fractions: count the favourable outcomes, divide by the total, apply Bayes' formula. This machinery works well enough for coins and dice, but it quietly falls apart the moment we ask about continuous quantities, infinite sample spaces, or the convergence of sequences of random variables. The reason is that probability is, at its core, a theory of measurement. The question "what is the probability of an event?" is really the question "what is the size of a set?", and answering it properly requires the same infrastructure that underpins the Lebesgue integral: -algebras, measures, and measurable functions.

This section rebuilds the foundations from that starting point. The payoff is not abstraction for its own sake. Once probability is grounded in measure theory, results like the law of large numbers and the central limit theorem become theorems about integrals, and the passage from discrete to continuous distributions requires no new definitions, only a change of reference measure. The Radon-Nikodym derivative, which we introduce at the end of this section, makes that change of reference precise and opens the door to the geometric perspective that runs through the rest of the post.

We begin with the two structures that any probabilistic model requires: a way to specify which events are observable, and a way to assign sizes to them.

Condition (iii) is the critical one: closure under countable unions, not merely finite ones. This is what separates -algebras from the simpler algebras of sets and is precisely the condition needed to take limits. Without it, we could not define convergence of events, and probability theory would be unable to make asymptotic statements.

Every set admits at least two sigma-algebras: the trivial sigma-algebra , in which only the certain and impossible events are observable, and the power set , in which every subset is observable. For uncountable spaces like , the power set is too large to support a useful measure, and we work instead with the Borel sigma-algebra , generated by the open intervals. The choice of sigma-algebra encodes what we can observe: a coarser sigma-algebra means less information.

The axioms are Kolmogorov's (1933), and their power lies in what they leave unspecified. They say nothing about equally likely outcomes, nothing about frequencies, nothing about subjective beliefs. They say only that probability is a normalised, countably additive set function. Every interpretation of probability (frequentist, Bayesian, decision-theoretic) is a commitment about what means; the measure-theoretic framework is agnostic.

With the probability space in hand, we can formalise what it means to assign a numerical value to each outcome.

The measurability condition [eq:eq:measurability] ensures that asking "what is the probability that falls in ?" is a well-posed question: the preimage belongs to , so can evaluate it. A function that is not measurable can produce sets whose probability is undefined, which is precisely the pathology that arises when one tries to assign probabilities to all subsets of (the Vitali construction).

The three-stage construction mirrors the Lebesgue integral exactly, because that is what expectation is. The only difference from an analyst's perspective is that the underlying measure is a probability measure rather than the Lebesgue measure. This identification is the central simplification of the measure-theoretic approach: every theorem about Lebesgue integrals (dominated convergence, monotone convergence, Fubini-Tonelli) applies immediately to expectations.

The following theorem illustrates why countable additivity, rather than mere finite additivity, is the right axiom. It lets us compute probabilities of limits from limits of probabilities.

Let us see all of these definitions at work in the simplest non-trivial probability space.

The Bernoulli space is minimal, but it already contains the essential structure: a sample space, a -algebra encoding what we can observe, a measure encoding our beliefs, and a random variable converting outcomes to numbers. Every probability model, no matter how complex, is built from the same ingredients.

1.1. Changing the reference: the Radon-Nikodym derivative

So far we have worked with a single probability measure . In practice we constantly compare two measures on the same measurable space: a prior and a posterior, a model and the data-generating distribution, a risk-neutral and a physical measure. The Radon-Nikodym derivative makes such comparisons precise.

1.2. The geometric thread: families as curves

Consider the family of all Bernoulli distributions parametrised by . Each value of defines a probability measure on , and these measures live in the probability simplex

The Bernoulli family traces a curve in this simplex: . The simplex is the ambient space of all distributions on a finite sample space, and parametric families are curves or surfaces within it. This geometric perspective, viewing families of distributions as geometric objects, will deepen throughout this post. For now, note that the curve passes through every point in , so the Bernoulli family is a full parametrisation of the simplex. Richer sample spaces will give higher-dimensional simplices, and the parametric families we study will typically be lower-dimensional submanifolds.

A natural question arises: which point on this curve is "most uncertain"? The answer leads to entropy.

from sympy import symbols, log, simplify, Rational
p = symbols('p', positive=True)
# Expectation of Ber(p)
EX = 0*(1-p) + 1*p
print(f"E[X] = {EX}")
# Variance
VarX = simplify(EX*(1-EX))  # E[X^2] - (E[X])^2 = p - p^2
print(f"Var(X) = {VarX}")
# Shannon entropy (base 2)
H = simplify(-p*log(p, 2) - (1-p)*log(1-p, 2))
print(f"H(p) = {H}")
# Value at p = 1/2
print(f"H(1/2) = {H.subs(p, Rational(1,2))} bits")

2. Estimation theory

With the measure-theoretic foundations in place, we turn to the central question of statistics: given observed data, what can we infer about the mechanism that generated it? This requires a formal notion of a statistical model and a systematic theory of estimation.

2.1. Statistical models and likelihood

The starting point is to specify what we mean by a "family of possible data-generating processes." A statistical model formalises this idea by indexing probability measures with a parameter.

The parametric assumption, that the family is smoothly indexed by a finite-dimensional parameter, is what makes classical estimation theory tractable. The dominating measure is typically Lebesgue measure (for continuous models) or counting measure (for discrete models). Its role is to allow us to work with densities rather than abstract measures.

Given a parametric model, an observation singles out a function on the parameter space: how "compatible" is each with what we saw?

The likelihood is emphatically not a probability distribution over : it does not integrate to one, and it carries no probabilistic interpretation on the parameter space (that would require a prior, which we defer to the Bayesian section). It is simply a score of compatibility. The log-likelihood is preferred in practice because products become sums, which are numerically stable and analytically convenient.

The most natural estimation principle is to choose the parameter value that makes the observed data most likely.

The MLE is a function of the data, hence a random variable. Its properties (consistency, efficiency, asymptotic normality) are among the deepest results in estimation theory, and we will establish them below after developing the necessary machinery.

2.2. Sufficient statistics

Not all functions of the data are equally useful. Some compress the data without losing any information about the parameter; others discard relevant information. The notion of sufficiency makes this precise.

Verifying sufficiency directly from the conditional distribution is often impractical. The following factorisation criterion provides a clean, checkable condition.

2.3. Exponential families

Sufficiency becomes particularly elegant in the context of exponential families, which encompass most of the standard parametric models and admit a rich geometric structure.

The factorisation theorem immediately tells us that is sufficient for . The log-partition function is convex (since it is the logarithm of a moment-generating function), and its derivatives encode the moments of the sufficient statistic.

The following result shows that the log-partition function is a generating function for the cumulants of the sufficient statistic.

2.4. Fisher information

The curvature of the log-likelihood at the true parameter value quantifies how much information the data carries about . This is formalised by the Fisher information.

from sympy import symbols, log, diff, simplify, Function
p, x = symbols('p x')
# Log-likelihood for single Bernoulli observation
logp = x*log(p) + (1-x)*log(1-p)
score = diff(logp, p)
print(f"Score: {score}")
# Fisher info: -E[d^2/dp^2 log p] = -E[second derivative]
d2 = diff(logp, p, 2)
# E[x] = p, E[1-x] = 1-p
I_p = simplify(-d2.subs(x, p))  # E[d2] by replacing x with E[x]=p
print(f"I(p) = {I_p}")

from sympy import symbols, log, diff, simplify, pi, sqrt, Matrix
x, mu, sigma = symbols('x mu sigma', real=True, positive=True)
# Log-density of N(mu, sigma^2)
logp = -log(sigma*sqrt(2*pi)) - (x - mu)**2 / (2*sigma**2)
# Score components
s_mu = diff(logp, mu)
s_sig = diff(logp, sigma)
print(f"Score (mu): {simplify(s_mu)}")
print(f"Score (sigma): {simplify(s_sig)}")
# Fisher info: I_ij = E[s_i * s_j]
# Use E[x] = mu, E[(x-mu)^2] = sigma^2, E[(x-mu)^4] = 3*sigma^4
from sympy import Rational
I_11 = simplify(s_mu**2)
I_11 = I_11.subs((x-mu)**2, sigma**2)  # take expectation
I_22 = simplify(s_sig**2)
I_22 = I_22.subs((x-mu)**4, 3*sigma**4).subs((x-mu)**2, sigma**2)
print(f"I_11 = E[s_mu^2] = {simplify(I_11)}")
print(f"I_22 = E[s_sigma^2] = {simplify(I_22)}")

The Fisher information admits a compelling geometric interpretation that will become central in Section 4.

2.5. Bias, variance, and the Cramer-Rao bound

Before establishing the optimality of the MLE, we need a framework for comparing estimators. The fundamental decomposition of estimation error separates systematic error (bias) from random fluctuation (variance).

The MSE decomposition reveals a fundamental tension: reducing bias (e.g. by shrinkage) may increase variance, and vice versa. The "best" estimator depends on which trade-off one is willing to accept. Within the class of unbiased estimators, however, there is a sharp lower bound on the variance.

The Cramer-Rao bound is tight in exponential families: for any exponential family, the MLE of achieves the bound exactly. Outside exponential families, the bound may not be attainable by any unbiased estimator, but the MLE still achieves it asymptotically.

2.6. Asymptotic efficiency of the MLE

The deepest result in classical estimation theory is that the MLE is, in a precise sense, the best estimator one can construct from large samples.

This theorem closes the circle: the Fisher information sets a fundamental limit on estimation precision (Cramer-Rao), and the MLE saturates that limit for large samples. The connection to geometry is immediate. The asymptotic covariance of the MLE is the inverse of the Fisher information matrix, which we will identify in Section 4 as the inverse of the Riemannian metric on the statistical manifold. Estimation theory and differential geometry are, in this sense, two perspectives on the same structure.

3. Inference and testing

Estimation tells us what to believe about a parameter; testing tells us what we can rule out. Hypothesis testing provides a formal framework for decisions under uncertainty, and its mathematical structure flows directly from the likelihood theory we have developed.

The central question shifts from "what is the best estimate of ?" to "is a particular value of compatible with the data?" To answer this, we need a decision procedure with quantifiable error guarantees.

3.1. The testing framework

Any decision procedure that divides the sample space into "reject" and "fail to reject" regions can make two kinds of mistakes. The entire theory of testing is built around controlling these errors.

A good test has small (rarely rejects a true null) and large power (reliably detects true alternatives). These goals are in tension: shrinking the critical region reduces but also reduces power. The Neyman-Pearson framework resolves this by fixing and maximising power.

Rather than choosing a significance level in advance, we can summarise the evidence against with a single number.

The uniformity statement is worth emphasising: if is true, the p-value is no more likely to fall below 0.05 than below any other threshold of the same width. This is what makes the p-value a valid calibration of evidence.

3.2. Optimal tests

Given the trade-off between Type I and Type II errors, a natural question arises: among all tests of a given size, which has the greatest power? The Neyman-Pearson lemma gives a complete answer for simple hypotheses.

The lemma tells us that the likelihood ratio is the sufficient statistic for distinguishing between two simple hypotheses. All the information relevant to the decision is captured in .

The Neyman-Pearson lemma handles simple vs simple hypotheses. For composite alternatives (e.g., ), we need a stronger condition on the statistical model.

Exponential families have monotone likelihood ratio in their natural sufficient statistics, so one-sided UMP tests exist for all standard exponential family models.

3.3. Large-sample testing

When exact UMP tests are unavailable (as with two-sided or multi-parameter problems), the likelihood ratio test provides an asymptotically optimal alternative.

Wilks' theorem is remarkably practical: it provides an asymptotic null distribution that does not depend on the true parameter value, only on the dimensions of the full and restricted parameter spaces. This universality is why the likelihood ratio test is the default in most applied settings.

3.4. Confidence regions

Testing and estimation are two sides of the same coin. The connection runs through confidence regions.

The construction of confidence regions often relies on quantities whose distributions are parameter-free.

The fundamental connection between testing and confidence regions is captured by the following duality.

This duality means that every test at level generates a confidence region (by inverting the acceptance region), and every confidence region generates a family of tests (by checking membership). The two frameworks carry exactly the same information.

3.5. Examples

3.6. Computational verification

The coverage probability of a confidence interval is a frequentist guarantee: over repeated samples, the interval should contain the true parameter at least of the time. We can verify this directly.

import numpy as np
np.random.seed(42)
p_true, n, alpha = 0.6, 30, 0.05
z = 1.96
n_sim = 10000
covers = 0
for _ in range(n_sim):
  x = np.random.binomial(n, p_true)
  p_hat = x / n
  se = np.sqrt(p_hat * (1 - p_hat) / n)
  lo, hi = p_hat - z * se, p_hat + z * se
  if lo <= p_true <= hi:
      covers += 1
print(f"Wald CI coverage: {covers/n_sim:.4f} (nominal: {1-alpha:.2f})")
print(f"(n={n}, p={p_true}, {n_sim} simulations)")

The Wald interval is known to undercover for moderate and extreme , a phenomenon that motivates alternatives like the Wilson interval. The gap between nominal and actual coverage is itself a statistical quantity that we can study.

The power curve of a test visualises how detection probability varies with the true parameter value. At , the power equals the size ; away from , it should increase toward 1.

import numpy as np
from scipy import stats
n = 20
alpha = 0.05
# Two-sided exact binomial test of H0: p = 0.5
p_values = np.linspace(0.01, 0.99, 50)
power = np.zeros_like(p_values)
for i, p in enumerate(p_values):
  # Power = P(reject H0 | true p)
  pw = 0
  for k in range(n + 1):
      pval = stats.binom_test(k, n, 0.5, alternative='two-sided')
      if pval <= alpha:
          pw += stats.binom.pmf(k, n, p)
  power[i] = pw
print("p      power")
for p, pw in zip(p_values[::5], power[::5]):
  print(f"{p:.2f}   {pw:.4f}")
print(f"\nPower at p=0.5: {power[24]:.4f} (should be ~alpha={alpha})")
print(f"Power at p=0.8: {power[39]:.4f}")

The power curve is symmetric about (since the test is two-sided) and increases monotonically as moves away from the null. The rate of increase depends on : larger samples yield steeper power curves, concentrating the "uncertain" region into a narrower band around .

4. The geometric synthesis

Throughout this post, we have spoken informally of the "curvature" of the log-likelihood, the "geometry" of the parameter space, and the "shape" of confidence regions. It is time to make these metaphors precise. The parameter space of a statistical model carries a natural Riemannian structure, the Fisher metric, that encodes the intrinsic geometry of statistical inference. In this section, we construct this geometry from first principles.

4.1. The manifold structure

4.2. The Riemannian structure

A manifold on its own has no notion of distance or angle. That structure is imposed by a metric.

Where does such a metric come from in statistics? The answer is the Fisher information.

We must verify that this definition actually yields a Riemannian metric: the key property is positive-definiteness.

The next theorem is the reason the Fisher metric deserves to be called geometric: it does not depend on the coordinate system.

4.3. Connection, geodesics, and curvature

A Riemannian metric determines a canonical way to differentiate vector fields along curves: the Levi-Civita connection. Its components in local coordinates are the Christoffel symbols.

4.4. The Fisher metric as infinitesimal KL divergence

The deepest justification for the Fisher metric is not algebraic but information-theoretic: it arises as the leading-order term in the Taylor expansion of the Kullback-Leibler divergence.

4.5. Worked examples: the geometry of familiar models

The Bernoulli manifold

The Normal manifold

We now compute the full Riemannian geometry of this surface: Christoffel symbols, geodesics, and curvature.

Christoffel symbols. The inverse metric is , , . Applying the Christoffel formula with coordinates , we note that all metric components depend only on (not on ), so any derivative . The nonzero derivatives are:

Computing each Christoffel symbol (all unlisted ones vanish):

from sympy import symbols, Matrix, simplify, Rational, diff
mu, sigma = symbols('mu sigma', positive=True)
# Fisher metric for N(mu, sigma): g = diag(1/sigma^2, 2/sigma^2)
g = Matrix([[1/sigma**2, 0], [0, 2/sigma**2]])
g_inv = g.inv()
coords = [mu, sigma]
print("Fisher metric g:")
print(g)
print("\nChristoffel symbols:")
for k in range(2):
  for i in range(2):
      for j in range(i, 2):
          val = 0
          for l in range(2):
              val += Rational(1,2) * g_inv[k,l] * (
                  diff(g[j,l], coords[i]) +
                  diff(g[i,l], coords[j]) -
                  diff(g[i,j], coords[l])
              )
          val = simplify(val)
          if val != 0:
              print(f"  Gamma^{k+1}_{{{i+1}{j+1}}} = {val}")

Gaussian curvature. We compute using the definition of the Riemann tensor:

Wait: let us be precise with the index convention. Using , we want with :

(All other -terms vanish because the relevant Christoffel symbols are zero.) Therefore:

Lowering the index: . The determinant of the metric is . The Gaussian curvature is:

The Normal statistical manifold has constant negative Gaussian curvature . This is the hallmark of hyperbolic geometry: triangles have angle deficit, parallel geodesics diverge, and circles grow exponentially.

from sympy import symbols, Matrix, simplify, Rational, diff
mu, sigma = symbols('mu sigma', positive=True)
g = Matrix([[1/sigma**2, 0], [0, 2/sigma**2]])
g_inv = g.inv()
coords = [mu, sigma]
# Compute Gamma^k_ij
def christoffel(k, i, j):
  val = 0
  for l in range(2):
      val += Rational(1,2) * g_inv[k,l] * (
          diff(g[j,l], coords[i]) + diff(g[i,l], coords[j]) - diff(g[i,j], coords[l]))
  return simplify(val)
# Riemann tensor R^1_212
R1_212 = (diff(christoffel(0,1,1), coords[0])
       - diff(christoffel(0,0,1), coords[1]))
for m in range(2):
  R1_212 += christoffel(0,0,m)*christoffel(m,1,1) - christoffel(0,1,m)*christoffel(m,0,1)
R1_212 = simplify(R1_212)
# R_1212 = g_1l R^l_212
R_1212 = simplify(g[0,0]*R1_212)
detg = simplify(g.det())
K = simplify(R_1212 / detg)
print(f"R^1_212 = {R1_212}")
print(f"R_1212 = {R_1212}")
print(f"det(g) = {detg}")
print(f"Gaussian curvature K = {K}")

Geodesics. The geodesic equations for the Normal manifold are, writing and :

These are the geodesic equations of the Poincare half-plane (with the coordinate rescaling reducing to the standard form). The solutions are of two types: vertical lines (corresponding to pure changes in variance) and semicircles centred on the -axis ( boundary). In statistical terms, the geodesic between two Normal distributions curves upward through higher-variance regions, because at larger the distributions are more spread out and therefore closer together in Fisher distance.

4.6. The Bayesian-frequentist synthesis

The geometric viewpoint does more than provide elegant computations: it reveals that the Bayesian and frequentist frameworks are two coordinate systems on the same Riemannian manifold. Their convergence at large sample sizes is a geometric theorem.

This theorem is the bridge. At large , Bayesian credible sets and frequentist confidence sets coincide, and both are determined by the Fisher geometry.

The visualisation above shows the Bernstein-von Mises theorem in action. At small , the prior matters: the Bayesian credible interval (gold) and frequentist confidence interval (blue) diverge. At large , both converge to the same interval, centred at the MLE, with width determined by the Fisher information. The geometric viewpoint explains why: both frameworks respect the Fisher metric. The frequentist approach conditions on and averages over data; the Bayesian approach conditions on data and integrates over . But the Fisher metric governs the precision of both. Where they differ is in what they condition on, not in the geometry they inhabit.

4.7. Looking ahead

The Fisher metric is the beginning, not the end, of information geometry. The full theory, developed by Amari and others, introduces -connections: a one-parameter family of connections on the statistical manifold, of which the Levi-Civita connection is the member. The and connections (the e-connection and m-connection) are dual with respect to the Fisher metric, and exponential families are dually flat with respect to this pair. This dual structure governs not just distance but divergence, projection, and the geometry of model selection. That is the subject of the sequel.

References

1. Kolmogorov, A. N. (1933). Grundbegriffe der Wahrscheinlichkeitsrechnung. Springer. English translation: Foundations of the Theory of Probability, 2nd ed., Chelsea, 1956. DOI
2. Fisher, R. A. (1922). On the Mathematical Foundations of Theoretical Statistics. Philosophical Transactions of the Royal Society A, 222(594-604), 309-368. DOI
3. Fisher, R. A. (1925). Theory of Statistical Estimation. Mathematical Proceedings of the Cambridge Philosophical Society, 22(5), 700-725. DOI
4. Neyman, J. and Pearson, E. S. (1933). On the Problem of the Most Efficient Tests of Statistical Hypotheses. Philosophical Transactions of the Royal Society A, 231(694-706), 289-337. DOI
5. Rao, C. R. (1945). Information and the Accuracy Attainable in the Estimation of Statistical Parameters. Bulletin of the Calcutta Mathematical Society, 37, 81-91. Link
6. Cramer, H. (1946). Mathematical Methods of Statistics. Princeton University Press. Publisher
7. Amari, S. and Nagaoka, H. (2000). Methods of Information Geometry. Translations of Mathematical Monographs, Vol. 191, AMS. DOI
8. Amari, S. (2016). Information Geometry and Its Applications. Applied Mathematical Sciences, Vol. 194, Springer. DOI
9. Lehmann, E. L. and Casella, G. (1998). Theory of Point Estimation, 2nd ed. Springer. DOI
10. Lehmann, E. L. and Romano, J. P. (2005). Testing Statistical Hypotheses, 3rd ed. Springer. DOI
11. van der Vaart, A. W. (1998). Asymptotic Statistics. Cambridge University Press. DOI
12. Jeffreys, H. (1946). An Invariant Form for the Prior Probability in Estimation Problems. Proceedings of the Royal Society A, 186(1007), 453-461. DOI
13. Le Cam, L. (1986). Asymptotic Methods in Statistical Decision Theory. Springer. DOI
14. do Carmo, M. P. (1992). Riemannian Geometry. Birkhauser. DOI
15. Efron, B. (1975). Defining the Curvature of a Statistical Problem (with Applications to Second Order Efficiency). Annals of Statistics, 3(6), 1189-1242. DOI
16. Kullback, S. and Leibler, R. A. (1951). On Information and Sufficiency. Annals of Mathematical Statistics, 22(1), 79-86. DOI