Random Samples, Joint Distributions, Sampling Statistics, and Moment Generating Techniques

At the core of statistical inference lies the mathematical formulation of a random sample. In empirical science, data collection generates a sequence of observed numerical values resulting from some physical, biological, or social mechanism. In the language of probability theory, these observations are modeled as random variables defined on an underlying probability space $(\Omega, \mathcal{F}, P)$, where $\Omega$ is the sample space, $\mathcal{F}$ is the associated $\sigma$-algebra of events, and $P$ is a probability measure.

Let $X_1, X_2, \dots, X_n$ be random variables mapping the sample space $\Omega$ into the real line $\mathbb{R}$. We state that the collection of random variables $\mathbf{X} = (X_1, X_2, \dots, X_n)$ constitutes a random sample of size $n$ from a population characterized by the cumulative distribution function $F(x; \boldsymbol{\theta})$ if two fundamental probabilistic conditions are simultaneously satisfied:

1. Identical Distribution: Each individual random variable $X_i$ possesses the exact same marginal probability distribution governed by the distribution function $F(x; \boldsymbol{\theta})$:

2. Mutual Stochastic Independence: The random variables $X_1, X_2, \dots, X_n$ are mutually independent, meaning that for any arbitrary sequence of Borel measurable sets $B_1, B_2, \dots, B_n \in \mathcal{B}(\mathbb{R})$, the joint probability factorizes as the product of the marginal probabilities:

This paired property is universally denoted by the acronym i.i.d. (independent and identically distributed). The statistical consequence of mutual independence is profound: the realization of any sub-collection of observations imparts zero probabilistic information regarding the behavior or values of the remaining observations, conditional on the true parameter $\boldsymbol{\theta}$.

To understand the scope of independence, it is illuminating to contrast it with the weaker condition of exchangeability. A sequence of random variables $X_1, \dots, X_n$ is said to be exchangeable if for every permutation $\pi$ of the indices $\{1, 2, \dots, n\}$, the joint probability distribution of the permuted vector $(X_{\pi(1)}, \dots, X_{\pi(n)})$ is identical to that of the original vector $(X_1, \dots, X_n)$. While every i.i.d. sequence is necessarily exchangeable due to the commutativity of real multiplication in product measures, exchangeable sequences need not be independent. A classic example is sampling without replacement from a finite population, where knowing that an earlier draw yielded a particular category slightly depletes that category, inducing negative covariance among draws while preserving permutation invariance. By de Finetti's Representation Theorem, any infinitely exchangeable sequence of random variables can be represented as a mixture of i.i.d. sequences with respect to some latent prior distribution, thereby anchoring the i.i.d. paradigm as the foundational mathematical building block of both frequentist and Bayesian statistical mechanics.

When the marginal distribution function $F(x; \boldsymbol{\theta})$ is dominated by a $\sigma$-finite reference measure $\mu$ on the measurable space $(\mathbb{R}, \mathcal{B}(\mathbb{R}))$ (such as the standard Lebesgue measure $\lambda$ in the continuous domain, or the discrete counting measure over a countable lattice), the Radon–Nikodym Theorem guarantees the existence of a marginal probability density function (or probability mass function) $f(x; \boldsymbol{\theta}) = \frac{dF}{d\mu}(x; \boldsymbol{\theta})$. Applying Tonelli's and Fubini's theorems to the product measure space, the joint distribution of the sample vector $\mathbf{X} = (X_1, \dots, X_n)$ is dominated by the product measure $\mu^n = \mu \times \dots \times \mu$ on $(\mathbb{R}^n, \mathcal{B}(\mathbb{R}^n))$, and its joint density function decomposes completely into the product:

This product decomposition serves as the mathematical engine for constructing likelihood functions, computing sampling moments, and establishing finite-sample as well as asymptotic properties of estimators.

Upon observing a realization $\mathbf{x} = (x_1, \dots, x_n)$ of the random sample $\mathbf{X}$, the raw data vector resides in an $n$-dimensional Euclidean space $\mathbb{R}^n$. In practical inference, an investigator rarely operates directly on the high-dimensional raw coordinates. Instead, the data vector is mapped into a lower-dimensional space through a function known as a statistic.

Formally, a statistic $T = T(\mathbf{X})$ is defined as a Borel measurable mapping:

where the functional form of $T$ contains no unknown parameters. The requirement that $T$ be free of unknown parameters distinguishes statistics from general random functions; for instance, the sample average $\frac{1}{n}\sum_{i=1}^n X_i$ is a legitimate statistic, whereas the quantity $\frac{1}{n}\sum_{i=1}^n (X_i - \mu)$ is a pivot or structural function, but not a statistic when the population mean $\mu$ is unknown.

From a measure-theoretic standpoint, the transformation $T$ induces a sub-$\sigma$-algebra $\sigma(T) \subseteq \mathcal{B}(\mathbb{R}^n)$, defined by the collection of inverse images:

The sub-$\sigma$-algebra $\sigma(T)$ represents the specific information concerning the underlying sample vector that can be discriminated by observing the value of $T$. Every statistic $T$ induces a partition of the sample space $\mathbb{R}^n$, where two sample vectors $\mathbf{x}$ and $\mathbf{y}$ belong to the same partition cell (equivalence class) if and only if $T(\mathbf{x}) = T(\mathbf{y})$. Data reduction consists precisely in collapsing each equivalence class into a single summary point, discarding within-cell variations while attempting to preserve between-cell variations that contain information about the unknown parameter $\boldsymbol{\theta}$.

Because the original sample vector $\mathbf{X}$ is a random vector governed by the probability measure $P_{\boldsymbol{\theta}}$, the statistic $T = T(\mathbf{X})$ is itself a random variable (or random vector if $m > 1$). The induced probability distribution of $T$, denoted $P_{\boldsymbol{\theta}}^T$ and defined on $(\mathbb{R}^m, \mathcal{B}(\mathbb{R}^m))$ by:

is termed the sampling distribution of the statistic $T$. The determination of exact sampling distributions is crucial for constructing exact confidence intervals, computing exact $p$-values in hypothesis testing, and evaluating risk functions in decision theory.

The most ubiquitous summary statistic for location is the sample mean, denoted $\bar{X}$ and defined as the arithmetic average of the sample coordinates:

The sample mean represents an affine projection of the $n$-dimensional sample vector onto the one-dimensional equicoordinate subspace spanned by the all-ones vector $\mathbf{1}_n = (1, 1, \dots, 1)^T \in \mathbb{R}^n$.

First Moment: Universality of Unbiasedness

Suppose that the population distribution possesses a finite first moment, meaning that $\mathbb{E}[|X_i|] < \infty$. Let $\mu = \mathbb{E}[X_i]$ denote the population mean. By the linearity of the mathematical expectation operator, the expectation of $\bar{X}$ is:

Notice that this derivation requires only the existence of the first moment $\mathbb{E}[|X_i|] < \infty$ and identical marginal expectations; it does not even require stochastic independence among the observations. Thus, the sample mean $\bar{X}$ is an unbiased estimator of the population mean $\mu$ across all probability distributions with finite mean.

Second Moment: Standard Error Decay under Independence

Suppose further that the population distribution possesses a finite second moment, so that $\text{Var}(X_i) = \sigma^2 < \infty$. Using the properties of the variance operator on sums of random variables:

Because the random variables $X_1, \dots, X_n$ are mutually independent, every pairwise covariance vanishes identically: $\text{Cov}(X_i, X_j) = 0$ for all $i \ne j$. Consequently, the double summation vanishes, leaving:

The positive square root of the sampling variance of $\bar{X}$ is the standard error of the mean:

This formula encapsulates the fundamental law of statistical precision: the variability of the sample mean diminishes inversely with the square root of the sample size, $\mathcal{O}(n^{-1/2})$. To halve the uncertainty surrounding the sample mean, an investigator must quadruple the sample size.

By Chebyshev's inequality, for any arbitrary positive tolerance $\epsilon > 0$:

Taking the limit as the sample size $n \to \infty$, the upper bound converges to zero:

This establishes that $\bar{X}$ converges in probability to the population mean $\mu$, proving the Weak Law of Large Numbers for any distribution with finite variance.

To quantify the dispersion or variability of the observed data around their empirical center, statistics relies on the sample variance, denoted $S^2$. The sample variance is defined with divisor $n - 1$:

The use of the denominator $n - 1$ instead of the sample size $n$ is known as Bessel's correction. To understand why this correction is necessary, consider the naive variance estimator using divisor $n$:

Let us derive the expectation of $\sum_{i=1}^n (X_i - \bar{X})^2$ by decomposing the deviation of each observation $X_i$ from the sample mean into two orthogonal components: the deviation of $X_i$ from the true population mean $\mu$, and the deviation of the sample mean $\bar{X}$ from $\mu$:

Squaring both sides of this identity yields:

Now, sum this equation across all observations $i = 1, 2, \dots, n$:

Notice that $\sum_{i=1}^n (X_i - \mu) = \sum_{i=1}^n X_i - n\mu = n\bar{X} - n\mu = n(\bar{X} - \mu)$. Substituting this into the cross-product term:

The third term is a constant summed $n$ times: $\sum_{i=1}^n (\bar{X} - \mu)^2 = n (\bar{X} - \mu)^2$. Combining these expressions yields the celebrated algebraic identity for the sum of squared deviations:

Now, take mathematical expectations across this identity:

By definition, $\mathbb{E}[(X_i - \mu)^2] = \text{Var}(X_i) = \sigma^2$ for every observation, so the first sum equals $n \sigma^2$. Furthermore, since $\bar{X}$ is an unbiased estimator of $\mu$, $\mathbb{E}[(\bar{X} - \mu)^2] = \text{Var}(\bar{X}) = \frac{\sigma^2}{n}$. Therefore:

This reveals that the naive estimator $S_n^2$ has expectation:

The naive estimator is systematically biased downwards; on average, it underestimates the true population variance by the exact amount $\frac{\sigma^2}{n}$. The geometric reason for this bias is intuitive: the deviations are measured from the sample mean $\bar{X}$, which is the unique center of gravity that minimizes the sum of squared distances $\sum_{i=1}^n (X_i - c)^2$ over all real choices $c$. Since $\bar{X}$ adapts to the specific sample, the sum of squared deviations from $\bar{X}$ is always strictly smaller than the sum of squared deviations from the true population mean $\mu$ (unless $\bar{X}$ happens to coincide exactly with $\mu$).

Dividing the sum of squared deviations by $n - 1$ instead of $n$ neutralizes this bias completely:

Thus, Bessel's correction guarantees that $S^2$ is an unbiased estimator of the population variance $\sigma^2$ for any arbitrary probability distribution possessing finite variance.

Linear combinations of random variables form the theoretical infrastructure of analysis of variance, linear regression, time series filtering, and multivariate analysis. Let $X_1, X_2, \dots, X_n$ be random variables with mean vector $\boldsymbol{\mu} = (\mu_1, \dots, \mu_n)^T$ and covariance matrix $\boldsymbol{\Sigma}$, where the $(i, j)$-th entry of $\boldsymbol{\Sigma}$ is $\sigma_{ij} = \text{Cov}(X_i, X_j)$.

Consider an arbitrary linear combination:

where $\mathbf{a} = (a_1, \dots, a_n)^T \in \mathbb{R}^n$ is a deterministic coefficient vector. By the linearity of the mathematical expectation:

To compute the variance of $Y$, we express the variance operator in quadratic form:

Since the inner product is a scalar, we can write:

Writing this quadratic form in scalar notation:

When the variables $X_1, \dots, X_n$ are mutually independent, the off-diagonal covariances vanish ($\sigma_{ij} = 0$ for all $i \ne j$), simplifying the quadratic form to a single diagonal sum:

Now, consider two distinct linear combinations:

The covariance between $Y_1$ and $Y_2$ is computed via bilinearity:

If the random variables $X_1, \dots, X_n$ are independent and identically distributed with common variance $\sigma^2$, then $\boldsymbol{\Sigma} = \sigma^2 \mathbf{I}_n$, where $\mathbf{I}_n$ is the $n \times n$ identity matrix. The covariance collapses to:

This provides a profound geometric connection between probability theory and Euclidean geometry: two linear combinations of i.i.d. variables are uncorrelated if and only if their coefficient vectors $\mathbf{a}$ and $\mathbf{b}$ are orthogonal in $\mathbb{R}^n$ with respect to the standard Euclidean inner product ($\mathbf{a}^T \mathbf{b} = 0$).

Orthogonal Decomposition of the Sample Space

To see this geometry in action, consider the decomposition of the sample vector $\mathbf{X} \in \mathbb{R}^n$ into the sample mean vector and the residual vector. Define the vector of sample means $\mathbf{\bar{X}} = \bar{X} \mathbf{1}_n = (\bar{X}, \bar{X}, \dots, \bar{X})^T$. The sample mean vector lies in the 1-dimensional subspace $L = \text{span}(\mathbf{1}_n)$. The residual vector is:

Notice that the sum of the components of $\mathbf{R}$ is identically zero: $\sum_{i=1}^n (X_i - \bar{X}) = 0$. In vector notation, this means:

Thus, the residual vector $\mathbf{R}$ is strictly orthogonal to the equicoordinate vector $\mathbf{1}_n$. It resides in the $(n - 1)$-dimensional orthogonal complement subspace $L^\perp$. The projection matrix that projects any vector $\mathbf{X}$ onto $L$ is:

The projection matrix that projects $\mathbf{X}$ onto $L^\perp$ is the centering matrix:

Both $\mathbf{P}$ and $\mathbf{M}$ are symmetric idempotent matrices: $\mathbf{P}^T = \mathbf{P}$, $\mathbf{P}^2 = \mathbf{P}$, $\mathbf{M}^T = \mathbf{M}$, and $\mathbf{M}^2 = \mathbf{M}$. Furthermore, they are mutually orthogonal: $\mathbf{P} \mathbf{M} = \mathbf{0}$. The sum of squared deviations can be written compactly as the squared Euclidean norm of the projected residual vector:

Because the rank of an idempotent matrix equals its trace, the dimension of the subspace $L^\perp$ where the residuals live is:

This geometric fact is the origin of the term "degrees of freedom": the $n$ original random variables are constrained to lie in an $(n - 1)$-dimensional hyperplane because their sum around the sample mean must equal zero.

While the direct calculation of sampling distributions via multivariate integrals can become intractable, transform methods provide an elegant algebraic machinery. Among these, the moment generating function (MGF) is the primary tool for analyzing sums and linear combinations.

Formal Definition and Convergence

Let $X$ be a random variable defined on $(\Omega, \mathcal{F}, P)$. The moment generating function of $X$, denoted $M_X(t)$, is defined as the mathematical expectation:

for all real numbers $t \in \mathbb{R}$ for which the expectation exists and is finite. We say that the moment generating function of $X$ exists if there exists a strictly positive constant $h > 0$ such that $M_X(t) < \infty$ for all $t \in (-h, h)$.

If $X$ is a continuous random variable with probability density function $f(x)$, the MGF is given by the integral:

If $X$ is a discrete random variable with probability mass function $p(x)$ supported on a countable set, the MGF is:

The condition that $M_X(t)$ is finite in an open interval around the origin $(-h, h)$ imposes a strong constraint on the decay rate of the distribution's tails. Specifically, for $t \in (0, h)$, the integral $\int_0^\infty e^{tx} f(x) \, dx$ converges only if the right tail of $f(x)$ decays faster than $e^{-tx}$ as $x \to +\infty$. Similarly, for $t \in (-h, 0)$, convergence requires the left tail to decay faster than $e^{-|t||x|}$ as $x \to -\infty$. Distributions with "heavy tails"—such as the Student's $t$, Pareto, and Cauchy distributions—fail this requirement, and their moment generating functions do not exist for any $t \ne 0$.

Differentiation and Moment Extraction

When $M_X(t)$ converges on an open interval $(-h, h)$, the function is infinitely differentiable (real analytic) on that domain. By the Dominated Convergence Theorem, we can interchange differentiation and expectation to any arbitrary order $k \in \mathbb{N}$:

Setting $t = 0$ collapses the exponential term to $e^0 = 1$, yielding the $k$-th raw moment $\mu_k' = \mathbb{E}[X^k]$:

Expanding $M_X(t)$ in a Maclaurin series reveals that the raw moments are directly encoded as the series coefficients:

The Uniqueness Theorem

The fundamental power of the MGF in distribution theory resides in its uniqueness.

Theorem (MGF Uniqueness Theorem): Let $X$ and $Y$ be two random variables with cumulative distribution functions $F_X$ and $F_Y$. If their moment generating functions $M_X(t)$ and $M_Y(t)$ both exist and satisfy:

then $F_X(u) = F_Y(u)$ at all continuity points $u \in \mathbb{R}$. That is, $X$ and $Y$ have identical probability distributions.

This theorem establishes an injective relationship between probability distributions and moment generating functions. If the MGF of an unknown statistic can be algebraically derived and shown to match the MGF of a known parametric distribution, the statistic must follow that exact distribution.

Algebraic Properties under Linear Operations

The MGF satisfies elementary algebraic properties under linear operations:

1. Affine Transformation: If $Y = aX + b$ for constants $a, b \in \mathbb{R}$, then:

2. Sums of Independent Random Variables: If $X_1, X_2, \dots, X_n$ are mutually independent random variables with respective MGFs $M_{X_1}(t), \dots, M_{X_n}(t)$, then the MGF of their sum $Y = \sum_{i=1}^n X_i$ is the product of their individual MGFs:

3. Sample Mean of i.i.d. Variates: For the sample mean $\bar{X} = \frac{1}{n} \sum_{i=1}^n X_i$ of an i.i.d. sample:

While the moment generating function provides a convenient tool when it exists, its failure to converge for heavy-tailed distributions limits its mathematical generality. To overcome this restriction, mathematical statistics deploys the characteristic function, which extends the real parameter $t$ into the complex plane by defining the argument along the imaginary axis.

The characteristic function of a random variable $X$ is defined for all real $t \in \mathbb{R}$ by:

where $i = \sqrt{-1}$ is the imaginary unit.

Because the integrand satisfies $|e^{itX}| = |\cos(tX) + i\sin(tX)| = \sqrt{\cos^2(tX) + \sin^2(tX)} = 1$, the integrand is universally bounded by 1 for all real numbers $t$ and $x$. Consequently:

Therefore, the characteristic function $\phi_X(t)$ exists for every probability distribution without exception, regardless of whether moments exist. The characteristic function possesses several universal structural properties:

1. $\phi_X(0) = \mathbb{E}[e^0] = 1$.
2. $\phi_X(-t) = \mathbb{E}[e^{-itX}] = \overline{\phi_X(t)}$, exhibiting Hermitian symmetry.
3. $\phi_X(t)$ is uniformly continuous on the entire real line $\mathbb{R}$.
4. If $Y = aX + b$, then $\phi_Y(t) = e^{itb} \phi_X(at)$.
5. If $X_1, \dots, X_n$ are independent, $\phi_{\sum X_i}(t) = \prod_{i=1}^n \phi_{X_i}(t)$.

Lévy's Inversion Formula

The characteristic function uniquely determines the underlying probability distribution, as established by Paul Lévy.

Theorem (Lévy's Inversion Theorem): Let $F(x)$ be the cumulative distribution function of $X$, and let $\phi(t)$ be its characteristic function. For any two real numbers $a < b$ that are continuity points of $F$:

Furthermore, if the characteristic function is absolutely integrable over the real line, $\int_{-\infty}^\infty |\phi(t)| \, dt < \infty$, then the distribution of $X$ is absolutely continuous with respect to Lebesgue measure, and its probability density function $f(x)$ is continuous, bounded, and given by the inverse Fourier transform:

This explicit inversion formula allows probability densities to be directly computed from characteristic functions, bridging distribution theory with harmonic analysis.

While moments provide intuitive polynomial expectations, their algebra under summation is complicated by binomial and multinomial coefficients. A vastly more natural algebraic representation is provided by cumulants, which are defined via the logarithm of the moment generating function.

The cumulant generating function (CGF) of a random variable $X$, denoted $K_X(t)$, is defined as:

assuming $M_X(t)$ exists in an open interval $(-h, h)$ around zero. The cumulants $\kappa_1, \kappa_2, \kappa_3, \dots$ are the coefficients of the Maclaurin series expansion of $K_X(t)$:

Equivalently, the $r$-th cumulant is given by the $r$-th derivative of the CGF evaluated at $t = 0$:

Mapping Between Cumulants and Central Moments

To derive the relationships between cumulants and moments, we expand $K_X(t) = \ln(1 + [M_X(t) - 1])$ using the Taylor series $\ln(1 + u) = u - \frac{u^2}{2} + \frac{u^3}{3} - \frac{u^4}{4} + \dots$, where $M_X(t) - 1 = \mu_1' t + \frac{\mu_2'}{2} t^2 + \frac{\mu_3'}{6} t^3 + \dots$:

1. First Cumulant ($\kappa_1$):

The first cumulant is the population mean.

2. Second Cumulant ($\kappa_2$):

The second cumulant is the population variance $\mu_2$.

3. Third Cumulant ($\kappa_3$):

Differentiating once more yields:

The third cumulant is the third central moment, capturing population skewness.

4. Fourth Cumulant ($\kappa_4$):

Computing the fourth derivative yields:

The fourth cumulant measures excess kurtosis. While the fourth central moment $\mu_4$ depends on variance, the fourth cumulant isolates the non-Gaussian tail behavior directly.

The Fundamental Property: Strict Additivity under Independence

The supreme advantage of cumulants is their behavior under independent summation. Let $X_1, \dots, X_n$ be mutually independent random variables, and let $Y = \sum_{i=1}^n X_i$. Taking the logarithm of the product of MGFs:

Equating coefficients of $\frac{t^r}{r!}$ on both sides demonstrates that for every integer order $r \ge 1$:

Cumulants are strictly additive under independent addition. There are no cross-terms or combinatorial expansions. To compute the $r$-th cumulant of a sum of independent random variables, one simply sums the individual $r$-th cumulants.

Marcinkiewicz's Theorem

For a Gaussian distribution $\mathcal{N}(\mu, \sigma^2)$, the MGF is $M(t) = \exp(\mu t + \frac{1}{2}\sigma^2 t^2)$, so the CGF is an exact quadratic polynomial:

Consequently, for any normal distribution, $\kappa_1 = \mu$, $\kappa_2 = \sigma^2$, and all higher-order cumulants vanish identically: $\kappa_r = 0$ for all $r \ge 3$.

A famous result by Józef Marcinkiewicz establishes that this is not an accident.

Theorem (Marcinkiewicz's Theorem): If the cumulant generating function $K(t)$ of a random variable is a polynomial in $t$, then the degree of that polynomial cannot exceed 2. That is, the polynomial must either be of degree 1 (a degenerate point-mass distribution) or of degree 2 (a Gaussian distribution).

Marcinkiewicz's Theorem reveals that no non-degenerate probability distribution can have its cumulants vanish beyond order 3, 4, or any higher finite order: either an infinite number of cumulants are non-zero, or only the first two are non-zero. This provides a deep mathematical characterization of the Gaussian distribution as the unique distribution having finite cumulant complexity.

When moment generating functions are inaccessible or when one requires the exact joint distribution of dependent variables, transformation of variables via Jacobian matrices and convolution integrals provides the fundamental constructive pathway.

The Multidimensional Change of Variables Theorem

Let $\mathbf{X} = (X_1, X_2)^T$ be a continuous bivariate random vector with joint probability density function $f_{\mathbf{X}}(x_1, x_2)$ supported on an open set $\mathcal{S} \subseteq \mathbb{R}^2$. Consider a transformation:

Assume that this mapping $\mathbf{g}: \mathcal{S} \to \mathcal{T} \subseteq \mathbb{R}^2$ is a bijective, continuously differentiable function with inverse mapping:

The Jacobian matrix of the inverse transformation is the $2 \times 2$ matrix of partial derivatives:

The Jacobian determinant, denoted $\det(J)$, captures the local area deformation factor under the coordinate change. By the multivariable change of variables theorem for Lebesgue integrals, the joint probability density function of $\mathbf{Y} = (Y_1, Y_2)^T$ on $\mathcal{T}$ is given by:

and zero elsewhere.

The Convolution Formula for Sums

To determine the probability density function of the sum $Y = X_1 + X_2$ of two continuous random variables, we introduce an auxiliary variable $Z = X_2$. The forward transformation is:

The inverse transformation is:

The Jacobian matrix of the inverse mapping is:

The determinant is $\det(J) = (1)(1) - (-1)(0) = 1$, so $|\det(J)| = 1$. The joint density of $(Y, Z)$ is therefore:

To obtain the marginal density of the sum $Y$, we integrate out the auxiliary variable $z$ over $\mathbb{R}$:

When $X_1$ and $X_2$ are independent, the joint density factorizes as $f_{X_1, X_2}(x_1, x_2) = f_{X_1}(x_1) f_{X_2}(x_2)$, yielding the classical convolution integral:

When both variables are non-negative ($X_1 \ge 0, X_2 \ge 0$), the integrand is non-zero only when $z \ge 0$ and $y - z \ge 0$, which restricts the integration limits to $z \in [0, y]$:

Quotient and Product Formulas

Similar Jacobian transformations yield standard formulas for ratios and products:

1. Ratio Transformation: Let $U = X_1 / X_2$ and $V = X_2$. The inverse is $X_1 = U V$ and $X_2 = V$. The Jacobian determinant is $|\det(J)| = |v|$. Integrating out $v$ yields the ratio density:

2. Product Transformation: Let $W = X_1 X_2$ and $V = X_2$. The inverse is $X_1 = W / V$ and $X_2 = V$. The Jacobian determinant is $|\det(J)| = 1/|v|$. Integrating out $v$ yields the product density:

A fundamental result connecting arbitrary continuous distributions with the standard uniform distribution is the Probability Integral Transform.

Theorem (Probability Integral Transform): Let $X$ be a continuous random variable whose cumulative distribution function $F_X(x) = P(X \le x)$ is strictly increasing and continuous on its support. Define the transformed random variable:

Then $U$ is distributed uniformly on the open unit interval $(0, 1)$:

Proof: Since $F_X$ is strictly increasing and continuous, its inverse function $F_X^{-1}: (0, 1) \to \mathbb{R}$ is strictly increasing and uniquely defined. For any $u \in (0, 1)$, the cumulative distribution function of $U$ is:

Applying the strictly increasing function $F_X^{-1}$ to both sides of the inequality preserves the inequality direction:

For $u \le 0$, $F_U(u) = 0$, and for $u \ge 1$, $F_U(u) = 1$. The CDF $F_U(u) = u$ on $(0, 1)$ is the exact CDF of a $\text{Uniform}(0, 1)$ random variable. Differentiating with respect to $u$ yields the constant density $f_U(u) = 1$ on $(0, 1)$. $\blacksquare$

The Inverse Transform Method for Simulation

The converse of the Probability Integral Transform provides the theoretical foundation for generating pseudo-random variates from any arbitrary target distribution.

Theorem (Inverse Transform Generation): Let $U \sim \text{Uniform}(0, 1)$, and let $F$ be any cumulative distribution function (continuous or discrete). Define the generalized inverse (quantile function):

Then the random variable $X = F^{-1}(U)$ has exact cumulative distribution function $F$.

This theorem guarantees that given a uniform random number generator, evaluating the inverse CDF $F^{-1}(U)$ generates exact, independent samples from any target distribution, anchoring computational Monte Carlo sampling in rigorous distribution theory.

In non-Gaussian parent populations, sample statistics frequently exhibit intricate statistical dependencies. While the sample mean $\bar{X}$ and sample variance $S^2$ are independent under normal sampling (as proved in Unit 3), they are generally correlated when sampling from skewed populations.

To derive the exact covariance between $\bar{X}$ and $S^2$, let $X_1, \dots, X_n$ be an i.i.d. random sample from a population with mean $\mu$, variance $\sigma^2$, and third central moment $\mu_3 = \mathbb{E}[(X - \mu)^3]$. We wish to evaluate:

Since $\mathbb{E}[\bar{X} - \mu] = 0$, the covariance simplifies to:

Recall the algebraic identity:

Multiplying by $\bar{X} - \mu = \frac{1}{n} \sum_{j=1}^n (X_j - \mu)$ yields:

Now, expand the first term:

When $j \ne i$, by independence:

When $j = i$, the expectation is the third central moment:

There are exactly $n$ such diagonal terms ($i = j \in \{1, \dots, n\}$). Therefore:

Next, we evaluate $\mathbb{E}[(\bar{X} - \mu)^3]$:

By mutual independence, any expectation where an index appears with power 1 vanishes because $\mathbb{E}[X_i - \mu] = 0$. The only terms that survive are those where all three indices coincide: $i = j = k$. There are $n$ such terms, each evaluating to $\mu_3$. Thus:

Substituting both expressions back into our covariance formula:

This is an exact, finite-sample result of fundamental importance:

We observe the following key implications:

1. If the population distribution is symmetric about its mean ($\mu_3 = 0$), such as in the Gaussian, Student's $t$, Uniform, or Laplace distributions, then $\text{Cov}(\bar{X}, S^2) = 0$. The sample mean and sample variance are completely uncorrelated.
2. If the population distribution is positively skewed ($\mu_3 > 0$), such as in the Exponential, Gamma, Lognormal, or Poisson distributions, then $\text{Cov}(\bar{X}, S^2) > 0$. When the sample mean is unusually large, the sample variance also tends to be unusually large.
3. If the population distribution is negatively skewed ($\mu_3 < 0$), such as in left-skewed Beta distributions, then $\text{Cov}(\bar{X}, S^2) < 0$.

This result underscores that zero correlation between $\bar{X}$ and $S^2$ is a universal feature of symmetric populations, but full stochastic independence between them is a unique signature reserved exclusively for Gaussian populations, as established by Geary's Theorem.

While the classical statistical framework presumes an infinite superpopulation or continuous data-generating mechanism yielding strictly independent observations, real-world survey methodology, ecological census tracking, and auditing procedures frequently sample from a finite, enumerated collection of $N$ identifiable population units. Let the population values be denoted by the finite sequence $\{y_1, y_2, \dots, y_N\}$, with fixed finite population mean $\mu_N = \frac{1}{N}\sum_{i=1}^N y_i$ and finite population variance $\sigma_N^2 = \frac{1}{N}\sum_{i=1}^N (y_i - \mu_N)^2$.

When drawing a sample of size $n < N$ without replacement (Simple Random Sampling Without Replacement, or SRSWOR), every possible subset of size $n$ has an identical probability of selection, equal to $1/\binom{N}{n}$. The resulting sample sequence $X_1, X_2, \dots, X_n$ is identically distributed: each $X_i$ has marginal probability distribution $P(X_i = y_k) = 1/N$ for all $k \in \{1, \dots, N\}$. Consequently, the marginal expectation and variance are identical across all draw indices:

However, because the sampling is executed without replacement, the observations are not mutually independent. Selecting a specific value $y_k$ on an earlier draw physically removes that value from the available urn, preventing its re-selection and thereby inducing a negative statistical dependency across all distinct coordinates. Specifically, for any pair $i \ne j$, the joint probability is:

with the probability vanishing when $k = m$. The joint expectation is:

Recognizing that $\sum_{k=1}^N y_k = N\mu_N$ and $\sum_{k=1}^N y_k^2 = N(\sigma_N^2 + \mu_N^2)$, the cross-product simplifies to:

Subtracting $\mathbb{E}[X_i]\mathbb{E}[X_j] = \mu_N^2$, we obtain the exact pairwise covariance between any two distinct draws:

This negative covariance reflects the physical conservation of units: observing an unusually large value depletes the remaining pool, making subsequent draws marginally smaller on average.

To determine the sampling variance of the sample mean $\bar{X} = \frac{1}{n}\sum_{i=1}^n X_i$ under SRSWOR, we sum the variances and covariances:

Factoring out $\frac{\sigma_N^2}{n}$:

Defining the population-adjusted variance $S_N^2 = \frac{N}{N-1}\sigma_N^2 = \frac{1}{N-1}\sum_{i=1}^N (y_i - \mu_N)^2$, this takes the classical form:

where $f = \frac{n}{N}$ is the sampling fraction, and the dimensionless factor:

is the Finite Population Correction (FPC). The FPC quantifies the exact reduction in sampling variability achieved when a substantial proportion of the entire population has been observed. When the population size $N \to \infty$ while the sample size $n$ remains fixed, the sampling fraction $f \to 0$, causing the FPC to converge to 1, gracefully recovering the standard infinite-population variance formula $\sigma^2 / n$. In the extreme limit where the entire population is sampled ($n = N$), the FPC collapses to zero, reflecting the deterministic fact that an exhaustive census possesses zero sampling error.

In experimental design, comparative clinical trials, and analysis of variance (ANOVA), researchers frequently investigate hypotheses concerning differences or weighted comparisons among group means. Such mathematical objects are termed linear contrasts.

Let $X_1, X_2, \dots, X_n$ be independent random variables with means $\mu_1, \dots, \mu_n$ and common variance $\sigma^2$. A linear combination:

is defined as a contrast if the scalar coefficients sum identically to zero:

Contrasts represent differences between subsets of variables. For instance, the simple pairwise difference $X_1 - X_2$ corresponds to $\mathbf{c} = (1, -1, 0, \dots, 0)^T$, and comparing the average of two treatment groups to a control group corresponds to $\mathbf{c} = (1/2, 1/2, -1, 0, \dots, 0)^T$.

The expectation of a contrast is the contrast of expectations:

If all group means are identical ($\mu_1 = \dots = \mu_n = \mu$), then $\mathbb{E}[C] = \mu \sum_{i=1}^n c_i = \mu(0) = 0$, regardless of the magnitude of $\mu$. The variance of the contrast is:

Now consider two distinct contrasts $C_1 = \mathbf{c}_1^T \mathbf{X}$ and $C_2 = \mathbf{c}_2^T \mathbf{X}$, where $\sum_{i=1}^n c_{1i} = 0$ and $\sum_{i=1}^n c_{2i} = 0$. The covariance between these two contrasts is:

Two contrasts are said to be orthogonal contrasts if their coefficient vectors are orthogonal in $\mathbb{R}^n$:

When contrasts are orthogonal, their covariance is zero. In a collection of $n$ observations, one can construct at most $n - 1$ mutually orthogonal contrasts that are simultaneously orthogonal to the equicoordinate vector $\mathbf{1}_n$. This geometric structure is operationalized through the Helmert Transformation.

The Helmert Transformation Matrix

The Helmert transformation is an explicit, constructive $n \times n$ orthogonal matrix $\mathbf{H}$ whose first row is proportional to $\mathbf{1}_n^T$, and whose remaining $n - 1$ rows constitute mutually orthogonal contrasts of normalized length 1:

By construction, $\mathbf{H}$ satisfies the orthogonality condition:

Let $\mathbf{Y} = \mathbf{H} \mathbf{X}$ denote the transformed random vector. Let us examine the individual components of $\mathbf{Y} = (Y_1, Y_2, \dots, Y_n)^T$:

1. The first component is:

The first transformed variable is an exact scalar multiple of the sample mean.

1. The remaining components $Y_2, Y_3, \dots, Y_n$ are mutually orthogonal contrasts that depend exclusively on pairwise and progressive differences among the observations, containing zero contribution from the grand sample mean.

Because $\mathbf{H}$ is an orthogonal transformation, it preserves the Euclidean inner product and Euclidean norm:

Now, expand the left side as the sum of squared transformed components:

Equating this to $\sum_{i=1}^n X_i^2$ yields the fundamental identity:

This establishes that the sum of squared deviations $(n - 1) S^2$ is identically equal to the sum of squares of the last $n - 1$ Helmert variables:

This algebraic decomposition provides the exact geometric foundation for proving the independence of the sample mean and sample variance in Gaussian populations, because the transformation maps the sample mean exclusively into $Y_1$, while mapping the sample variance exclusively into the remaining orthogonal coordinates $Y_2, \dots, Y_n$.

When sampling pairs of measurements $(X_{1i}, X_{2i})$ from a shared physical or biological subject (such as height and weight, parental and offspring trait values, or pre-test and post-test scores), the univariate framework generalizes to bivariate and multivariate distributions. The canonical baseline for joint continuous modeling is the bivariate normal distribution.

A bivariate random vector $\mathbf{X} = (X_1, X_2)^T$ follows a bivariate normal distribution, denoted $\mathbf{X} \sim \mathcal{N}_2(\boldsymbol{\mu}, \boldsymbol{\Sigma})$, with mean vector $\boldsymbol{\mu} = (\mu_1, \mu_2)^T$ and positive-definite covariance matrix:

where $\sigma_1 > 0$, $\sigma_2 > 0$, and the Pearson correlation coefficient $\rho \in (-1, 1)$ measures the degree of linear association between the two components. The determinant of the covariance matrix is:

The inverse covariance matrix (the precision matrix) is computed explicitly:

The quadratic form appearing in the exponent is:

Thus, the joint probability density function of $(X_1, X_2)$ is given for all $(x_1, x_2) \in \mathbb{R}^2$ by:

Conditional Distributions and Linear Regression

A celebrated structural property of the bivariate normal distribution is that its conditional distributions are strictly univariate normal. By completing the square in the exponent with respect to $x_2$, the joint density factorizes as $f_{X_1, X_2}(x_1, x_2) = f_{X_1}(x_1) f_{X_2 \mid X_1}(x_2 \mid x_1)$, where the marginal density is $X_1 \sim \mathcal{N}(\mu_1, \sigma_1^2)$ and the conditional density is:

This conditional law reveals two foundational mathematical properties:

1. Linearity of the Regression Function: The conditional expectation $\mathbb{E}[X_2 \mid X_1 = x_1] = \mu_2 + \beta (x_1 - \mu_1)$ is a strictly linear function of the conditioning value $x_1$, where the regression slope $\beta = \rho \frac{\sigma_2}{\sigma_1}$ represents the rate of change in $X_2$ per unit change in $X_1$. This provides the rigorous probabilistic justification for linear regression models.

2. Homoskedasticity: The conditional variance $\text{Var}(X_2 \mid X_1 = x_1) = \sigma_2^2 (1 - \rho^2)$ is completely constant; it does not depend in any way on the conditioning realization $x_1$. The variance of $X_2$ is reduced from its marginal value $\sigma_2^2$ by the factor $1 - \rho^2$, quantifying the exact proportion of variance explained by knowing $X_1$ ($R^2 = \rho^2$).

Sampling Distribution of the Sample Correlation Coefficient

In empirical investigations, the population correlation $\rho$ is unknown and must be estimated from a bivariate random sample $(X_{11}, X_{21}), \dots, (X_{1n}, X_{2n})$ via the Pearson sample correlation coefficient:

When the true population correlation is zero ($\rho = 0$), the exact sampling distribution of $R$ was established by Student (William Sealy Gosset). The transformed statistic:

follows an exact Student's $t$ distribution with $n - 2$ degrees of freedom:

However, when $\rho \ne 0$, the exact probability density function of $R$ (first derived by R. A. Fisher in 1915) is notoriously complicated, involving Gaussian hypergeometric functions ${}_2F_1$. Furthermore, the distribution of $R$ becomes highly skewed as $|\rho|$ approaches 1, converging to normality very slowly.

To resolve this inferential challenge, Fisher devised the famous Fisher $z$-transformation:

Fisher proved via Taylor expansion that the transformed statistic $Z$ stabilizes the asymptotic variance and converges rapidly to a Gaussian distribution:

Notice that the asymptotic variance $\frac{1}{n - 3}$ is completely independent of the unknown parameter $\rho$. This variance-stabilizing transformation allows analysts to construct robust, high-accuracy confidence intervals and hypothesis tests for correlation coefficients across arbitrary sample sizes.