Probability and Statistics

Probability is the mathematics of uncertainty, and statistics is the practice of learning about uncertain quantities from data. Computer vision needs both at every level. A pixel value is a noisy measurement of the light that reached the sensor; a detector’s output is a guess with some chance of being wrong; a classifier’s softmax output is read as a probability; a tracker carries a belief about where an object is; and every reported accuracy is an estimate from a finite test set. This article is an overview of the ideas these uses share: random variables and their distributions, joint and conditional probability, Bayes’ rule, expectation, variance and covariance, what happens when random quantities are added and averaged, and how quantities are estimated from samples. It links to the articles that develop each use, such as Gaussian Distribution and Image Noise, rather than repeating them.

Definition

A probability assigns to each event, a set of possible outcomes, a number between 0 and 1, such that a certain event has probability 1 and the probabilities of mutually exclusive events add. A random variable is a quantity whose value is uncertain, such as the count recorded by a pixel in the next exposure or the class of the object in an image, and its distribution says how probability is spread over its possible values. Statistics runs the other way: given observed values (samples), it estimates properties of the distribution that produced them, such as its mean or covariance, and says how far those estimates can be trusted.

Intuition

Photograph a static, evenly lit grey card many times without moving the camera. A given pixel reads 121, then 118, then 124: it never repeats exactly, because the number of photons arriving during an exposure is itself random. The distribution of that pixel’s readings is what probability describes. Its centre (the expectation) is what the scene “really” looks like at that pixel, and its spread (the variance) is the noise. Collecting the readings and drawing a histogram gives an empirical distribution, which approaches the true one as more exposures are added; averaging the readings gives an estimate of the expectation that gets better the more frames are averaged.

Now point the camera at a person and ask whether a given pixel shows skin. The pixel’s colour is evidence, but some background objects, such as wooden furniture, have similar colours. Knowing how skin and background colours are distributed, and how much of a typical image is skin, Bayes’ rule turns the observed colour into the probability that the pixel is skin. Most of vision’s use of probability is one of these two moves: describing the variability of measurements, and reasoning backwards from a measurement to its likely cause.

Formal Definition

The notation follows the site’s shared convention: $p(x)$ for a probability mass or density function, $P(\cdot)$ for the probability of an event, and $\mathbb{E}$, $\operatorname{Var}$ and $\operatorname{Cov}$ for expectation, variance and covariance. Vectors are bold and are columns.

Random variables and distributions

A discrete random variable $X$ takes values in a finite or countable set: a pixel’s 8-bit value, a photon count, a class label. Its probability mass function (PMF) gives the probability of each value, $p(x) = P(X = x)$, with $0 \le p(x) \le 1$ and $\sum_x p(x) = 1$ (Murphy, 2022, §2.2.1). The normalized image histogram of an image is an empirical PMF of its intensities: the fraction of pixels at each level.

A continuous random variable takes real values, such as the irradiance at a point or a sub-pixel keypoint position. The probability of any single exact value is zero, so the distribution is described by the cumulative distribution function (CDF) and its derivative, the probability density function (PDF):

$$F(x) = P(X \le x), \qquad p(x) = \frac{dF(x)}{dx}, \qquad P(a < X \le b) = \int_a^b p(x)\, dx = F(b) - F(a).$$

A density satisfies $p(x) \ge 0$ and $\int p(x)\, dx = 1$ (Bishop, 2006, §1.2.1; Murphy, 2022, §2.2.2). Nothing bounds it above: the probability of a short interval of width $\delta x$ is about $p(x)\,\delta x$, so a distribution concentrated in a narrow range has densities larger than 1. The CDF is defined for discrete variables too, where it is the cumulative histogram. Its inverse gives quantiles: the median is the value at which $F$ reaches $0.5$ (Murphy, 2022, §2.2.2.3).

Joint, marginal and conditional distributions

Two random variables, such as the colour $\mathbf{c}$ of a pixel and its class $Y$ (skin or background), have a joint distribution $p(\mathbf{c}, y)$. Summing or integrating out one variable gives the marginal of the other, the sum rule; dividing the joint by a marginal gives a conditional distribution, from which follows the product rule (Bishop, 2006, §1.2):

$$p(\mathbf{c}) = \sum_{y} p(\mathbf{c}, y), \qquad p(y \mid \mathbf{c}) = \frac{p(\mathbf{c}, y)}{p(\mathbf{c})}, \qquad p(\mathbf{c}, y) = p(\mathbf{c} \mid y)\, p(y).$$

$p(\mathbf{c} \mid y = \text{skin})$ is the distribution of colours among skin pixels; $p(y = \text{skin} \mid \mathbf{c})$ is the fraction of pixels of colour $\mathbf{c}$ that are skin. They are different quantities, and confusing them is one of the most common errors in reasoning with probabilities.

$X$ and $Y$ are independent when the joint factorizes, $p(x, y) = p(x)\, p(y)$, so that observing one says nothing about the other: $p(y \mid x) = p(y)$. They are conditionally independent given $Z$ when $p(x, y \mid z) = p(x \mid z)\, p(y \mid z)$ (Murphy, 2022, §2.2.4). Vision models lean heavily on both assumptions. “White” noise means the noise at different pixels is uncorrelated, and models usually go further and take it to be independent, and many models treat pixels or features as conditionally independent given a class or a scene, which makes their joint probability a product of simple factors.

Bayes’ rule

Writing the joint both ways with the product rule and dividing gives Bayes’ rule (Bishop, 2006, §1.2; Murphy, 2022, §2.3):

$$p(y \mid \mathbf{c}) = \frac{p(\mathbf{c} \mid y)\, p(y)}{p(\mathbf{c})} = \frac{p(\mathbf{c} \mid y)\, p(y)}{\sum_{y'} p(\mathbf{c} \mid y')\, p(y')}.$$

$p(y)$ is the prior, what is known about the class before seeing the pixel; $p(\mathbf{c} \mid y)$, viewed as a function of $y$ for the observed colour, is the likelihood; and the result is the posterior. The denominator only normalizes, so in words, posterior $\propto$ likelihood $\times$ prior. The rule is named after Thomas Bayes, whose essay on the problem of inverse probability was found among his papers after his death and communicated to the Royal Society by Price, in a letter read on 23 December 1763; it is printed in volume 53 of the Philosophical Transactions, dated 1763 (Bayes, 1763), which Bishop (2006, §1.2.3) gives as published in 1764.

For two classes, the posterior exceeds $\tfrac{1}{2}$ exactly when the likelihood ratio exceeds the inverse ratio of the priors:

$$p(\text{skin} \mid \mathbf{c}) > \tfrac{1}{2} \iff \frac{p(\mathbf{c} \mid \text{skin})}{p(\mathbf{c} \mid \text{background})} > \frac{p(\text{background})}{p(\text{skin})}.$$

Assigning every input to the class with the largest posterior minimizes the probability of a misclassification (Bishop, 2006, §1.5.1). When skin covers 7% of an image, a colour must be about 13 times more likely under skin than under background before the minimum-error decision calls it skin. The same structure appears whenever vision reasons from evidence to causes. Deciding whether a large image gradient is a real edge or noise weighs how often noise alone produces such a gradient against how often edges do, and how common edges are; recovering 3D shape from a 2D image, a problem with many solutions consistent with the same image, uses a prior to prefer plausible shapes (Murphy, 2022, §2.3.3). Treating unknown model parameters themselves as random variables with priors, and updating them with Bayes’ rule, is Bayesian inference; the Bayes filter applies it recursively to a moving object’s state.

Expectation, variance, covariance and correlation

The expectation of a function $f$ of a random variable is its probability-weighted average (Bishop, 2006, §1.2.2):

$$\mathbb{E}[f(X)] = \sum_x f(x)\, p(x) \quad \text{or} \quad \int f(x)\, p(x)\, dx .$$

The mean $\mu = \mathbb{E}[X]$ and the variance $\operatorname{Var}[X] = \mathbb{E}[(X - \mu)^2] = \mathbb{E}[X^2] - \mu^2$ summarize location and spread; the standard deviation $\sigma$ is the square root of the variance and has the same units as $X$. Expectation is linear, $\mathbb{E}[aX + b] = a\mathbb{E}[X] + b$, while $\operatorname{Var}[aX + b] = a^2 \operatorname{Var}[X]$ (Murphy, 2022, §2.2.5).

For two variables, the covariance measures how they vary together, and the correlation coefficient normalizes it to lie in $[-1, 1]$ (Murphy, 2022, §3.1):

$$\operatorname{Cov}[X, Y] = \mathbb{E}\big[(X - \mathbb{E}[X])(Y - \mathbb{E}[Y])\big], \qquad \rho = \frac{\operatorname{Cov}[X, Y]}{\sqrt{\operatorname{Var}[X] \operatorname{Var}[Y]}} .$$

Correlation measures the degree of linear relationship: $\rho = \pm 1$ only when $Y$ is an exact linear function of $X$, and independent variables have zero covariance, but zero covariance does not imply independence (Murphy, 2022, §3.1.2–3.1.3).

For a random vector $\mathbf{x} \in \mathbb{R}^d$, such as a pixel’s RGB colour or a feature descriptor, the covariance matrix collects the variances on its diagonal and the covariances off it:

$$\Sigma = \operatorname{Cov}[\mathbf{x}] = \mathbb{E}\big[(\mathbf{x} - \boldsymbol{\mu})(\mathbf{x} - \boldsymbol{\mu})^\top\big], \qquad \operatorname{Cov}[A\mathbf{x} + \mathbf{b}] = A \Sigma A^\top .$$

$\Sigma$ is symmetric and positive semidefinite (Murphy, 2022, §3.1.1). The transformation rule propagates uncertainty through any linear map, such as a colour-space conversion or the motion model of a Kalman filter, and the eigenvectors of $\Sigma$ (see linear algebra) give the directions of greatest variation, the starting point of principal component analysis.

Properties

Common distributions in vision

A handful of distributions cover most of what vision models need:

• Bernoulli and binomial. A Bernoulli variable is 1 with probability $\theta$ and 0 otherwise: a pixel is foreground or not, a detection is correct or not. The number of successes in $n$ independent Bernoulli trials is binomial, with mean $n\theta$ and variance $n\theta(1-\theta)$ (Bishop, 2006, §2.1). The number of correctly classified images in a test set is modelled this way, which is where the standard error of accuracy comes from.
• Categorical. One of $K$ classes, with probabilities $\theta_1, \dots, \theta_K$ summing to one (Murphy, 2022, §2.5). A classifier’s softmax output is a categorical distribution over classes, and the cross-entropy loss is the negative log-probability it assigns to the true label.
• Poisson. A count $k \in \{0, 1, 2, \dots\}$ with $p(k) = e^{-\lambda} \lambda^k / k!$ (Murphy, 2022, §12.2.3). It models the number of photons a pixel collects during an exposure; its variance equals its mean $\lambda$, so photon shot noise grows with brightness (Image Noise develops the sensor model).
• Uniform. Every value in an interval is equally likely. Rounding a measurement to integer levels produces an error that is close to uniform over half a level either side, which is how quantization noise is modelled.
• Gaussian. The bell-shaped density defined by a mean and a covariance, the default model for measurement noise and for uncertainty in estimation and tracking, and the subject of Gaussian Distribution, which covers its geometry, the Mahalanobis distance and its closure properties.

Sums and averages of random variables

For any random variables, the mean of a sum is the sum of the means. For variances, the covariances enter (Murphy, 2022, §2.8.4):

$$\operatorname{Var}[X_1 + X_2] = \operatorname{Var}[X_1] + \operatorname{Var}[X_2] + 2 \operatorname{Cov}[X_1, X_2],$$

so for independent variables, variances add. The full distribution of a sum of independent variables is the convolution of their distributions (Murphy, 2022, §2.8.5).

Averaging $n$ independent measurements with mean $\mu$ and variance $\sigma^2$ therefore gives, by the two rules above,

$$\mathbb{E}[\bar{x}] = \mu, \qquad \operatorname{Var}[\bar{x}] = \frac{\sigma^2}{n}, \qquad \operatorname{std}[\bar{x}] = \frac{\sigma}{\sqrt{n}} .$$

The variance of the average shrinks to zero as $n$ grows, so the probability that the sample average differs from the expectation by more than any fixed amount goes to zero: the law of large numbers, which Grinstead and Snell (2006, §8.1) prove in exactly this way for independent trials with finite variance. It is why an expectation can be approximated by the average of samples drawn from the distribution, with the approximation becoming exact as $n \to \infty$ (Bishop, 2006, §1.2.2). Averaging $n$ frames of a static scene reduces the noise standard deviation by $\sqrt{n}$; four times as many frames halve the noise. The central limit theorem adds that, for independent terms with a common distribution of finite variance, the distribution of the suitably normalized sum or average approaches a Gaussian as $n$ grows, whatever that common distribution is (Murphy, 2022, §2.8.6; Grinstead and Snell, 2006, §9.2–9.3); versions with weaker conditions allow terms that are not identically distributed (Grinstead and Snell, 2006, §9.2; Bishop, 2006, §2.3). This is one reason Gaussian noise models work and why averaged measurements are treated as Gaussian.

Sampling and estimation

Statistics treats a data set as samples from an unknown distribution and computes estimators: functions of the samples that estimate a property of the distribution. With $n$ samples stored one per row in a matrix $X \in \mathbb{R}^{n \times d}$, as scikit-learn stores them, the usual estimators of the mean and covariance are

$$\hat{\boldsymbol{\mu}} = \frac{1}{n} \sum_{i=1}^{n} \mathbf{x}_i, \qquad \hat{\Sigma} = \frac{1}{n - 1} \sum_{i=1}^{n} (\mathbf{x}_i - \hat{\boldsymbol{\mu}})(\mathbf{x}_i - \hat{\boldsymbol{\mu}})^\top .$$

An estimator is itself a random variable: drawn again, the data would give a different estimate. Its bias is the difference between its expectation and the true value, and its variance measures how much it changes from one data set to another; the expected squared error decomposes as (Murphy, 2022, §4.7.6)

$$\mathbb{E}\big[(\hat{\theta} - \theta)^2\big] = \operatorname{Var}[\hat{\theta}] + \operatorname{bias}(\hat{\theta})^2 .$$

The sample mean is unbiased. The sample variance computed with denominator $n$ is biased low, with expectation $\frac{n-1}{n}\sigma^2$, because the deviations are measured from the sample mean, which sits closer to the samples than the true mean does; dividing by $n - 1$ removes the bias (Bishop, 2006, §1.2.4; Murphy, 2022, §4.7.6.1). The denominator matters for small samples, such as a variance computed from a 3 × 3 patch, and the libraries disagree on the default: NumPy’s np.var and np.std divide by $n$ (ddof=0), while np.cov divides by $n - 1$ (NumPy 2.3.5 documentation). np.cov also treats each row as a variable by default, so a data matrix with one sample per row needs rowvar=False.

The most widely used general recipe for building estimators is maximum likelihood estimation: choose the parameter values under which the observed data are most probable. For Gaussian data it gives the sample mean and the variance with denominator $n$ (Bishop, 2006, §1.2.4). Adding a prior over the parameters and maximizing the posterior instead gives maximum a posteriori (MAP) estimation.

How much an estimate can be trusted is expressed by its standard error, the standard deviation of its sampling distribution, and by confidence intervals built from it. Statistical Significance in Evaluation applies these tools to the scores of vision models, where the test set is the sample.

Examples

Photon counts and frame averaging

The code renders the expected photon counts of a dim scene, a disc at 8 photons per pixel on a background at 2, draws Poisson counts at every pixel, compares the background’s empirical PMF with the Poisson PMF, and then averages $N$ independent frames.

import numpy as np
from scipy import stats

rng = np.random.default_rng(0)

# Expected photon counts of a dim synthetic scene: a disc (8 photons per pixel)
# on a background (2 photons per pixel), 128 x 128 pixels.
yy, xx = np.mgrid[:128, :128]
disc = (yy - 64) ** 2 + (xx - 64) ** 2 < 30 ** 2
lam = np.where(disc, 8.0, 2.0)
bg = ~disc                                    # background pixels, lambda = 2

# One frame: an independent Poisson count at every pixel.
frame = rng.poisson(lam)
k = np.arange(6)
print("background counts k:", k)
print("empirical PMF:      ", np.round(np.bincount(frame[bg])[:6] / bg.sum(), 3))
print("Poisson(2) PMF:     ", np.round(stats.poisson.pmf(k, 2.0), 3))

# Average N independent frames; measure the noise over background pixels.
print(f"\n{'N':>3} {'mean':>6} {'std':>6} {'sqrt(2/N)':>9} {'skewness':>8}")
for n in [1, 4, 16, 64]:
    avg = rng.poisson(lam, size=(n, *lam.shape)).mean(axis=0)
    r = avg[bg]
    print(f"{n:>3} {r.mean():>6.3f} {r.std():>6.3f} {np.sqrt(2 / n):>9.3f} "
          f"{stats.skew(r):>8.3f}")

Output:

background counts k: [0 1 2 3 4 5]
empirical PMF:       [0.136 0.272 0.27  0.184 0.086 0.034]
Poisson(2) PMF:      [0.135 0.271 0.271 0.18  0.09  0.036]

  N   mean    std sqrt(2/N) skewness
  1  1.978  1.402     1.414    0.722
  4  2.000  0.709     0.707    0.315
 16  2.005  0.354     0.354    0.172
 64  1.998  0.175     0.177    0.087

The histogram of the 13,575 background pixels from a single frame already matches the Poisson PMF to within half a percentage point at every count: across pixels that share the same expected value, the image itself is a sample from the distribution. Averaging leaves the mean at 2, since the average is unbiased, while the standard deviation follows $\sqrt{\lambda / N}$, halving each time the number of frames quadruples. The skewness, the standardized third moment, which is zero for a symmetric distribution such as the Gaussian, also roughly halves with each quadrupling: a single dim frame has lopsided noise, with counts piled up near zero and a long tail of bright pixels, and the average of many frames looks increasingly Gaussian, as the central limit theorem predicts. Averaging works only for noise that changes between frames; a static scene and camera are assumed, and fixed-pattern noise, which repeats in every frame, is not reduced.

The covariance of pixel colours

The code builds an RGB image with an orange square on a blue-grey background. A random shading factor per pixel brightens or darkens all three channels together, and independent noise is added to each channel. It computes the covariance matrix of the pixel colours with denominator $n - 1$, checks it against np.cov, and compares the correlation of red and blue over the whole image with that inside each region.

import numpy as np

rng = np.random.default_rng(0)

# A synthetic 96 x 128 RGB image: an orange square on a blue-grey background.
# Each region has a base colour, a shading factor that brightens or darkens
# all three channels together, and independent sensor noise per channel.
h, w = 96, 128
obj = np.zeros((h, w), bool)
obj[24:72, 40:88] = True
base = np.where(obj[..., None], [200.0, 120.0, 60.0], [70.0, 90.0, 110.0])
shade = 1 + 0.15 * rng.standard_normal((h, w, 1))
img = base * shade + rng.normal(0, 5, (h, w, 3))

X = img.reshape(-1, 3)                  # n x 3 data matrix, one pixel per row
n = X.shape[0]
Xc = X - X.mean(axis=0)
S = Xc.T @ Xc / (n - 1)                 # sample covariance, denominator n - 1
print("n =", n)
print("covariance (n - 1):\n", np.round(S, 1))
print("matches np.cov(X, rowvar=False):", np.allclose(S, np.cov(X, rowvar=False)))
print("np.cov(X[:10]).shape (rows taken as variables):", np.cov(X[:10]).shape)
print("correlation:\n", np.round(np.corrcoef(X, rowvar=False), 2))

# Correlation of red and blue within each region versus the whole image.
r = [np.corrcoef(X[m], rowvar=False)[0, 2] for m in (obj.ravel(), ~obj.ravel())]
print(f"correlation R-B: object {r[0]:.2f}, background {r[1]:.2f}")

Output:

n = 12288
covariance (n - 1):
 [[2886.2  819.9 -802.1]
 [ 819.9  373.6  -19.5]
 [-802.1  -19.5  639.5]]
matches np.cov(X, rowvar=False): True
np.cov(X[:10]).shape (rows taken as variables): (10, 10)
correlation:
 [[ 1.    0.79 -0.59]
 [ 0.79  1.   -0.04]
 [-0.59 -0.04  1.  ]]
correlation R-B: object 0.86, background 0.86

The hand-computed matrix agrees with np.cov once rowvar=False is set; without it, NumPy treats the pixels as the variables and returns a pixel-by-pixel matrix. With 12,288 pixels, the choice between $n$ and $n - 1$ changes the entries by less than 0.01%.

Within each region, red and blue are strongly positively correlated, because shading scales all channels together. Over the whole image their correlation is negative, because the image mixes two populations: the object has high red and low blue, the background the opposite, and the difference between the two region means dominates the shading within each. This is the law of total variance, which splits a variance into the average variance within groups plus the variance of the group means (Murphy, 2022, §2.2.5.4), at work on covariances. A covariance computed over a whole image, or over a feature set pooled from several classes, can describe the mixture rather than any of its parts.

Estimating noise from small patches: bias and variance

Noise level estimators often work on small flat patches. The code cuts a flat noisy image (noise standard deviation 5, variance 25) into 10,000 patches of 3 × 3 pixels, estimates the variance in each patch with both denominators, and measures each estimator’s bias, spread and mean squared error.

import numpy as np

rng = np.random.default_rng(0)

# A flat grey 300 x 300 image with Gaussian noise of standard deviation 5,
# cut into 10,000 patches of 3 x 3 pixels. Each patch estimates the noise
# variance (true value 25) from its own 9 pixels.
img = 128 + rng.normal(0, 5, (300, 300))
patches = img.reshape(100, 3, 100, 3).transpose(0, 2, 1, 3).reshape(-1, 9)

print(f"std of patch means: {patches.mean(axis=1).std():.3f}  (5 / 3 = {5 / 3:.3f})")
for ddof in [0, 1]:
    v = patches.var(axis=1, ddof=ddof)  # divide by 9 - ddof
    print(f"ddof={ddof}: mean {v.mean():6.2f}  bias {v.mean() - 25:6.2f}  "
          f"std {v.std():5.2f}  MSE {np.mean((v - 25) ** 2):6.1f}")

Output:

std of patch means: 1.676  (5 / 3 = 1.667)
ddof=0: mean  22.23  bias  -2.77  std 11.11  MSE  131.1
ddof=1: mean  25.01  bias   0.01  std 12.50  MSE  156.2

The patch means scatter with the standard deviation $\sigma / \sqrt{9}$ that the averaging rule predicts. Dividing by $n = 9$ underestimates the variance by the factor $8/9$, which is what the bias formula gives ($25 \times 8/9 = 22.2$); dividing by $n - 1 = 8$ removes the bias. Unbiased is not the same as accurate, though: each 9-pixel estimate is so noisy that its spread is half the true value, and the biased estimator, whose estimates are all scaled down by $8/9$, has the smaller mean squared error. Which estimator is better depends on what the error costs, a trade-off between bias and variance that runs through statistics and machine learning. A noise estimate that must be reliable pools many patches rather than trusting any one.

Bayes’ rule for skin-colour pixels

The code renders a scene with skin-coloured discs on a background of a blue wall, a grey floor and a light-brown wooden panel whose colours overlap those of skin. It learns class-conditional colour histograms $p(\mathbf{c} \mid \text{class})$ over 512 colour bins (8 levels per channel) and the prior $P(\text{skin})$ from a labelled training image, then classifies the pixels of a new image of the same kind with Bayes’ rule, and compares the result with the rule that ignores the prior and compares likelihoods alone. One pseudo-count is added to every bin, so that a colour never seen in training does not get probability zero. Five random scenes are run.

import numpy as np

def render(rng, h=160, w=160):
    """Synthetic scene: five skin-coloured discs (label 1) on a background of
    a blue wall, a grey floor and a light-brown wooden panel (label 0)."""
    yy, xx = np.mgrid[:h, :w]
    img = np.empty((h, w, 3))
    img[:] = [60, 80, 140]                                   # wall
    img[yy > 0.6 * h] = [110, 110, 105]                      # floor
    img[(xx > 0.65 * w) & (yy < 0.6 * h)] = [180, 125, 90]   # wooden panel
    label = np.zeros((h, w), bool)
    for _ in range(5):
        cy, cx = rng.uniform(20, h - 20, 2)
        label |= (yy - cy) ** 2 + (xx - cx) ** 2 < rng.uniform(8, 14) ** 2
    img[label] = [205, 145, 120]                             # skin
    img = img * (1 + 0.1 * rng.standard_normal((h, w, 1)))   # shading
    img += rng.normal(0, 8, img.shape)                       # sensor noise
    return np.clip(img, 0, 255), label

def colour_bin(img, bins=8):
    q = (img // (256 // bins)).astype(int)                   # 8 levels per channel
    return (q[..., 0] * bins + q[..., 1]) * bins + q[..., 2]  # one of 512 colours

def report(pred, truth):
    fp, fn = np.sum(pred & ~truth), np.sum(~pred & truth)
    return f"error {np.mean(pred != truth):.3f} (false skin {fp:4d}, missed {fn:4d})"

for seed in range(5):
    rng = np.random.default_rng(seed)
    train_img, train_lab = render(rng)
    test_img, test_lab = render(rng)

    # Class-conditional colour histograms p(c | class), with one pseudo-count
    # per bin so unseen colours keep a small probability, and the prior P(skin).
    b = colour_bin(train_img)
    counts = [np.bincount(b[train_lab == k], minlength=512) + 1.0 for k in (0, 1)]
    lik = [c / c.sum() for c in counts]
    prior = train_lab.mean()

    bt = colour_bin(test_img)
    joint1 = lik[1][bt] * prior
    post = joint1 / (joint1 + lik[0][bt] * (1 - prior))      # Bayes' rule
    print(f"seed {seed}: prior P(skin) {prior:.3f}, "
          f"{test_lab.sum()} skin pixels in the test image\n"
          f"  posterior > 0.5:      {report(post > 0.5, test_lab)}\n"
          f"  likelihood ratio > 1: {report(lik[1][bt] > lik[0][bt], test_lab)}")

Output:

seed 0: prior P(skin) 0.073, 2125 skin pixels in the test image
  posterior > 0.5:      error 0.057 (false skin   86, missed 1365)
  likelihood ratio > 1: error 0.072 (false skin 1717, missed  130)
seed 1: prior P(skin) 0.065, 1583 skin pixels in the test image
  posterior > 0.5:      error 0.043 (false skin   76, missed 1035)
  likelihood ratio > 1: error 0.067 (false skin 1623, missed   87)
seed 2: prior P(skin) 0.074, 2068 skin pixels in the test image
  posterior > 0.5:      error 0.052 (false skin   69, missed 1273)
  likelihood ratio > 1: error 0.063 (false skin 1483, missed  119)
seed 3: prior P(skin) 0.081, 2144 skin pixels in the test image
  posterior > 0.5:      error 0.056 (false skin   66, missed 1366)
  likelihood ratio > 1: error 0.062 (false skin 1486, missed  108)
seed 4: prior P(skin) 0.094, 1731 skin pixels in the test image
  posterior > 0.5:      error 0.046 (false skin  669, missed  512)
  likelihood ratio > 1: error 0.059 (false skin 1406, missed   92)

In all five scenes the posterior rule makes fewer errors than the likelihood-only rule: 4.3–5.7% of the 25,600 pixels against 5.9–7.2%. The two rules disagree only on colours whose likelihood ratio lies between 1 and the ratio of the priors, colours that both classes produce, which in these scenes come from skin and the wooden panel, whose base colours were chosen close together. Skin covers only 6.5–9.4% of each training image, so the posterior rule needs a likelihood ratio of about 10 to 14 before it calls a pixel skin and labels most of the shared colours background, while the likelihood-only rule calls them skin and marks 1,406–1,717 background pixels as skin. The minimum-error rule pays for its lower error by missing skin: 512 of 1,731 skin pixels with seed 4, and more than 60% of them with the other four seeds. Whether that is acceptable depends on what a missed skin pixel costs compared with a false one; moving the posterior threshold trades one error for the other, the trade-off that precision and recall and ROC curves describe. The classifier also inherits the histogram’s limits: it ignores where a pixel is and what its neighbours look like, and the prior learned from training images is only right if test images contain skin in similar proportions.

Common Misconceptions

• “$p(A \mid B)$ and $p(B \mid A)$ are about the same.” They can differ by orders of magnitude. Most skin pixels may be orange-brown, yet most orange-brown pixels in a scene full of wooden furniture are not skin; the prior links the two.
• “Uncorrelated means independent.” Zero correlation rules out only a linear relationship. A variable and its square can be uncorrelated while one determines the other (Murphy, 2022, §3.1.3).
• “A density is a probability, so it is at most 1.” A density can take any non-negative value; only its integral over a region is a probability. A Gaussian with a standard deviation of 0.1 pixels has a peak density of about 4 per pixel.
• “Averaging more frames removes all noise.” It reduces only the noise that is independent between frames, by the square root of the number of frames. Fixed-pattern noise, motion in the scene and systematic errors remain.
• “The covariance of the whole data set describes each part of it.” Pooled statistics mix the differences between groups with the variation within them, and can even reverse the sign of a correlation, as in the colour example.
• “A classifier’s softmax output is the probability that it is right.” It is a categorical distribution that the network was trained to produce. Whether its confidence matches how often it is actually right has to be checked on test images; expected calibration error measures the gap.

Where It Is Used

• Image formation and noise: Poisson shot noise, Gaussian read noise and quantization noise in image noise models; frame averaging and noise estimation; image denoising, which starts from an assumed noise distribution.
• Histograms and colour models: the image histogram as an empirical distribution, used for thresholding, colour-based detection and tracking appearance models.
• Estimation and geometry: least squares as the maximum-likelihood estimate under Gaussian noise; covariance matrices for the uncertainty of fitted parameters; the Gaussian distribution and Mahalanobis distance for gating in data association.
• Tracking: the Bayes filter and its Gaussian and sample-based forms, the Kalman filter and the particle filter.
• Learning: classifiers that output categorical distributions through the softmax, trained by minimizing the cross-entropy loss, which is maximum likelihood; Gaussian mixture models, principal component analysis and other models built on means and covariances; information theory, which supplies entropy and the KL divergence.
• Evaluation: binomial uncertainty of accuracy, confidence intervals and tests in statistical significance in evaluation, and calibration of predicted probabilities in expected calibration error.