Probability & Statistics 3: Statistical Estimation, Inequalities, and Limit Theorems

The first article in the series established probability spaces, conditional probability, and random variables; the second article in the series covered common distributions, numerical summaries, joint distributions, and correlation. Those two articles started from a probability model, assuming the distribution or parameters were already known; statistical inference faces the opposite problem: the researcher only has access to finite data, and wants to know the population distribution and unknown parameters that produced it.

This article uses one running set of concrete numbers throughout, distinguishing between a statistic, an estimator, and an estimate, introducing the method of moments and maximum likelihood estimation, and covering bias, variance, mean squared error, and consistency. Finally, it derives why the sample mean is stable starting from the Markov and Chebyshev inequalities, then explains exactly what question each of the law of large numbers, the central limit theorem, and standard error answers.

1. From the Population to a Random Sample

1.1 Population, Sample, and Observations

The population is the entire set of individuals or generating mechanism a research question is concerned with — for example, every part produced on a given production line, every request a service will ever receive in the future, or the complete distribution of a given random variable. The population mean $\mu$, population variance $\sigma^2$, and a Bernoulli success probability $p$ are all parameters; a parameter describes the population, and is usually unknown and fixed.

A researcher obtains $n$ observations from the population. Before observation, $X_1,\ldots,X_n$ are random variables; after observation, $x_1,\ldots,x_n$ are the concrete values already obtained. The distinction between uppercase and lowercase matters: $\bar X$ changes across repeated sampling and therefore has a sampling distribution, while $\bar x=3$ computed from one particular dataset is just a single number.

The most commonly used theoretical model is an independent and identically distributed (IID) random sample: every $X_i$ follows the same distribution, and different observations are mutually independent. IID isn’t a property data automatically has. Continuously measured network latency may have temporal dependence, multiple records generated by the same user may cluster together, and convenience sampling may fail to represent the population. When the sampling mechanism violates independence or representativeness, adding more data points doesn’t automatically eliminate a systematic bias.

1.2 Statistic, Estimator, and Estimate

A statistic is a function of the sample, and its formula must not contain any unknown parameter. The sample mean $\bar X=n^{-1}\sum_{i=1}^nX_i$, the sample maximum $\max_iX_i$, and the empirical distribution function are all statistics. $\bar X-\mu$ contains the unknown $\mu$, so it isn’t a statistic that can be computed directly from the data.

When a statistic is used to infer an unknown parameter, it’s called an estimator. For example, $\hat\mu=\bar X$ is an estimator of the population mean. The number obtained by plugging actual data into an estimator is called an estimate; if the five observations are $1,2,2,4,6$, then $\hat\mu=3$. The same estimation formula is a random variable before sampling, and produces a determinate result after sampling.

2. Computing Common Statistics from a Sample

2.1 Central Location, Proportion, Moments, and Variance

Consider five days’ worth of failure counts $x=(1,2,2,4,6)$. The sample mean is $\bar x=(1+2+2+4+6)/5=3$, and the sample median is 2. The median only depends on the middle position after sorting, whereas the mean is affected by every value; if the failure count 6 became 60, the median would still be 2, but the mean would increase from 3 to 13.8.

Binary data commonly use the sample proportion. If the inspection results for ten products are $0,1,0,0,1,0,0,0,1,0$, with 1 meaning defective, then the number of defects is $K=3$ and the sample proportion is $\hat p=K/n=0.3$. Since the mean of binary values also equals the proportion of 1s, a Bernoulli sample satisfies $\bar x=\hat p$ as well.

The $k$-th order sample raw moment is $m'_k=n^{-1}\sum_i x_i^k$. For the failure data, the first-order raw moment is $m'_1=3$, and the second-order raw moment is $m'_2=(1^2+2^2+2^2+4^2+6^2)/5=12.2$. The $k$-th order sample central moment is $m_k=n^{-1}\sum_i(x_i-\bar x)^k$; the second-order central moment is $m_2=[(-2)^2+(-1)^2+(-1)^2+1^2+3^2]/5=3.2$.

The $m_2$ with denominator $n$ describes the average squared deviation of the current five data points, and also shows up in the maximum likelihood estimator for the normal distribution. If the goal is an unbiased estimate of the population variance, $s^2=(n-1)^{-1}\sum_i(x_i-\bar x)^2=16/4=4$ is usually used instead. The denominators $n$ and $n-1$ correspond to different purposes, and shouldn’t be swapped out of habit. The sample standard deviation is $s=2$, in the same units as the original failure count.

2.2 Empirical Distribution, Quantiles, and Covariance

The empirical cumulative distribution function (empirical CDF) is defined as $F_n(t)=n^{-1}\sum_i\mathbf 1(x_i\le t)$ — that is, the proportion of the sample not exceeding $t$. For the failure data, $F_n(2)=3/5=0.6$ and $F_n(4)=4/5=0.8$. The empirical CDF preserves the full distributional information of a one-dimensional dataset without needing to first assume normal, Poisson, or some other parametric distribution.

Using the empirical quantile definition $q_p=\inf\{t:F_n(t)\ge p\}$, the failure data gives $q_{0.25}=2$, $q_{0.5}=2$, $q_{0.9}=6$. Different software may interpolate between adjacent order statistics differently, so quantile results for small samples can differ slightly. When reporting a quantile, note which definition or software default was used — don’t mistake a difference between interpolation rules for a computational error.

Paired data needs to describe how two variables vary together. Suppose the input for four trials is $x=(1,2,3,4)$ and the output is $y=(2,3,5,8)$, with means $\bar x=2.5$ and $\bar y=4.5$ respectively. The sample covariance is $s_{xy}=(n-1)^{-1}\sum_i(x_i-\bar x)(y_i-\bar y)=10/3\approx3.33$. A positive sign indicates that a larger input tends to be accompanied by a larger output, but the magnitude of covariance changes depending on the units of measurement.

Standardizing the covariance gives the sample correlation coefficient $r=s_{xy}/(s_xs_y)=10/\sqrt{105}\approx0.976$. A correlation coefficient close to 1 indicates a strong positive linear relationship in the data, but this alone doesn’t prove that the input causes the output to increase — common causes, the sampling method, and nonlinear relationships all need to be analyzed separately. The model-level definitions of joint distributions, marginal distributions, and correlation can be found back in the second article in the series.

3. Estimating Parameters with the Method of Moments and Maximum Likelihood

3.1 The Method of Moments Aligns Theoretical Moments with Sample Moments

The method of moments first writes out the relationship between a distribution’s moments and its parameters, then sets the theoretical moments equal to the corresponding sample moments. If $X\sim\operatorname{Bernoulli}(p)$, then $E[X]=p$; setting $E[X]=\bar X$ gives $\hat p_{\mathrm{MM}}=\bar X$. For the ten products above with three defects, the method-of-moments estimate is 0.3. The Poisson distribution satisfies $E[X]=\lambda$, so the Poisson rate parameter can likewise be estimated with $\hat\lambda_{\mathrm{MM}}=\bar X$.

If $X\sim N(\mu,\sigma^2)$, the first two raw moments satisfy $E[X]=\mu$ and $E[X^2]=\mu^2+\sigma^2$. Setting the theoretical moments equal to $m'_1$ and $m'_2$ respectively gives $\hat\mu_{\mathrm{MM}}=\bar X$ and $\hat\sigma^2_{\mathrm{MM}}=m'_2-\bar X^2=n^{-1}\sum_i(X_i-\bar X)^2$. The failure data yields $\hat\mu=3$ and $\hat\sigma^2=12.2-3^2=3.2$.

The method of moments is usually easy to compute, and also serves well as a starting point for numerical optimization, but it doesn’t necessarily make full use of the distribution’s information. Higher-order sample moments are also easily affected by extreme values. Which moments to choose, and whether the resulting estimate falls within a valid parameter range, both need to be checked against the specific model.

3.2 Maximum Likelihood Selects the Parameter That Best Explains the Data

Given observations $x_1,\ldots,x_n$, the likelihood function treats the joint PMF or PDF as a function of the parameter $\theta$, written $L(\theta;x)=\prod_i f(x_i;\theta)$. With the data held fixed, maximum likelihood estimation (MLE) searches for the parameter value that maximizes $L(\theta;x)$. Likelihood is not “the probability that the parameter is true”; in a frequentist model, the parameter isn’t a random variable, and $L(\theta_1;x)/L(\theta_2;x)$ only compares the relative support the same data provides for two parameter values.

For a Bernoulli sample, if there are $k$ successes out of $n$ trials, the likelihood is $L(p)=p^k(1-p)^{n-k}$. The product is usually rewritten as the log-likelihood $\ell(p)=k\log p+(n-k)\log(1-p)$, since taking the log doesn’t change the location of the maximum, and turns the product into an easier-to-handle sum. Setting the derivative $\ell'(p)=k/p-(n-k)/(1-p)$ equal to 0 gives $\hat p_{\mathrm{MLE}}=k/n$. The three defects out of ten products again gives $\hat p=0.3$.

For an IID normal sample, the log-likelihood is $\ell(\mu,\sigma^2)=-\frac n2\log(2\pi\sigma^2)-\frac{1}{2\sigma^2}\sum_i(x_i-\mu)^2$. Maximizing separately over $\mu$ and $\sigma^2$ gives $\hat\mu_{\mathrm{MLE}}=\bar x$ and $\hat\sigma^2_{\mathrm{MLE}}=n^{-1}\sum_i(x_i-\bar x)^2$. For the failure data, the normal MLE gives $\hat\mu=3$, $\hat\sigma^2=3.2$. The MLE for the normal variance uses denominator $n$, even though it’s the sample variance with denominator $n-1$ that’s unbiased — maximum likelihood and unbiasedness are two different criteria.

MLE’s answer depends on the complete model. If the observations are actually dependent on one another, a censoring mechanism is ignored, or count data is mistakenly treated as normal data, differentiating a formally correct expression can still produce a meaningless estimate. Before starting the optimization, you should clearly write out the distribution of each observation, the parameter range, and how the joint likelihood decomposes.

4. Judging Whether an Estimator Is Reliable

4.1 Bias, Variance, and Mean Squared Error

An estimator needs to be evaluated under “repeatedly drawing samples of the same size.” The bias of an estimate of $\theta$ is defined as $\operatorname{Bias}(\hat\theta)=E[\hat\theta]-\theta$; when the expectation equals the parameter, the estimator is called unbiased. For an IID sample, $E[\bar X]=\mu$, so the sample mean is unbiased; the Bernoulli sample proportion $\hat p=\bar X$ is likewise unbiased.

The sample variance with denominator $n$, $\tilde S^2=n^{-1}\sum_i(X_i-\bar X)^2$, satisfies $E[\tilde S^2]=(n-1)\sigma^2/n$, which systematically underestimates the population variance. The reason is that the same batch of data was already used to estimate $\bar X$ first, so the deviations are subject to one linear constraint, leaving only $n-1$ degrees of freedom. Multiplying by $n/(n-1)$ gives $S^2=(n-1)^{-1}\sum_i(X_i-\bar X)^2$, which satisfies $E[S^2]=\sigma^2$.

Unbiasedness doesn’t mean an estimate is necessarily close to the true value. An estimator’s mean squared error is $\operatorname{MSE}(\hat\theta)=E[(\hat\theta-\theta)^2]=\operatorname{Var}(\hat\theta)+\operatorname{Bias}(\hat\theta)^2$. MSE penalizes both sampling fluctuation and systematic bias at the same time, so an estimator with very small bias but very large variance can end up worse than a stable estimator carrying a small amount of bias. Small-sample Bernoulli problems sometimes use Laplace smoothing, $\tilde p=(K+1)/(n+2)$; this estimator is usually biased, but it avoids directly producing an estimate of 0 or 1 when $K=0$ or $K=n$, and can lower the MSE outside extreme parameter regions.

4.2 Consistency and Standard Error

If $\hat\theta_n$ converges in probability to $\theta$ as the sample size increases, the estimator is called consistent — that is, $P(|\hat\theta_n-\theta|>\varepsilon)\to0$ for every $\varepsilon>0$. Consistency describes a large-sample limit, and doesn’t guarantee accuracy at any particular finite sample size; unbiasedness describes the sampling average at each fixed sample size, and doesn’t guarantee small variance either. The two properties answer different questions.

The standard deviation of an estimator’s sampling distribution is called the standard error (SE). If an IID sample has population variance $\sigma^2$, then $\operatorname{Var}(\bar X)=\sigma^2/n$, so $\operatorname{SE}(\bar X)=\sigma/\sqrt n$; since $\sigma$ is usually unknown, it’s typically replaced with $s$, giving the estimated standard error $s/\sqrt n$. The sample standard deviation $s$ describes the spread of individual observations, while the standard error $s/\sqrt n$ describes the sampling fluctuation of the sample mean — the two values can’t be used interchangeably.

5. From Probability Inequalities to Limit Theorems

5.1 The Markov and Chebyshev Inequalities

If a random variable $Y\ge0$ and $E[Y]<\infty$, the Markov inequality gives $P(Y\ge a)\le E[Y]/a$, for $a>0$. If request latency is non-negative with a mean of 100 ms, the Markov inequality can only guarantee $P(Y\ge500\text{ ms})\le0.2$. The bound can be quite loose, but the conclusion doesn’t require knowing the full distribution of the latency.

Applying the Markov inequality to the non-negative variable $(X-\mu)^2$ gives the Chebyshev inequality: $P(|X-\mu|\ge a)\le\sigma^2/a^2$. Substituting a multiple of the standard deviation, $a=k\sigma$, turns this into $P(|X-\mu|\ge k\sigma)\le1/k^2$. For instance, any distribution with finite variance satisfies $P(|X-\mu|\ge2\sigma)\le1/4$; the Chebyshev inequality assumes nothing about normality, so it can’t be casually swapped out for the normal distribution’s empirical 95% rule.

5.2 The Weak Law of Large Numbers, the Strong Law of Large Numbers, and the Central Limit Theorem

For an IID sample with mean $\mu$ and variance $\sigma^2<\infty$, the sample mean satisfies $E[\bar X_n]=\mu$ and $\operatorname{Var}(\bar X_n)=\sigma^2/n$. Applying the Chebyshev inequality to $\bar X_n$ gives $P(|\bar X_n-\mu|\ge\varepsilon)\le\sigma^2/(n\varepsilon^2)$. The right-hand side approaches 0 as $n$ increases, so $\bar X_n$ converges in probability to $\mu$. This derivation gives the weak law of large numbers under the condition of finite variance.

The strong law of large numbers states a stronger property about the sample path: under suitable conditions, $P(\lim_{n\to\infty}\bar X_n=\mu)=1$. The weak law says that, at a fixed large sample size, the probability of deviating from $\mu$ approaches 0; the strong law says that, except for a set of sample paths with probability 0, the running average computed along the same infinite sampling sequence converges. Almost-sure convergence doesn’t mean every conceivable sequence converges, and it doesn’t mean any finite-sample average is exactly equal to $\mu$ either.

The law of large numbers answers whether an estimator approaches its target, but doesn’t describe the shape of the remaining error. If an IID sample has finite mean $\mu$ and positive finite variance $\sigma^2$, the classical central limit theorem (CLT) gives $\sqrt n(\bar X_n-\mu)/\sigma\Rightarrow N(0,1)$. So, for large samples, this can be approximately written as $\bar X_n\sim N(\mu,\sigma^2/n)$, with the typical error scale being $1/\sqrt n$. Quadrupling the sample size roughly halves the standard error; multiplying the sample size by a hundred shrinks the standard error to roughly a tenth of its original value.

Take a Bernoulli sample of $n=400$ as an example: if $\hat p=0.30$ is observed, plugging this into the estimated standard error gives $\widehat{\operatorname{SE}}(\hat p)=\sqrt{0.3\times0.7/400}\approx0.0229$. Under the normal approximation, the approximate 95% margin of error is $1.96\times0.0229\approx0.0449$, corresponding to an interval of roughly $[0.255,0.345]$. This interval is an approximate result based on the sampling distribution; extremely small samples, very rare success events, strongly dependent data, or infinite-variance distributions can all cause the normal approximation to break down.

The central limit theorem also doesn’t mean the original data’s distribution itself becomes normal. What approaches normality is the centered and rescaled sample sum or sample mean. The law of large numbers is concerned with whether the location of $\bar X_n$ is stable, while the central limit theorem is concerned with the scale and approximate shape of the fluctuation of $\bar X_n-\mu$. To extend both limit theorems to partial sums, random walks, and Brownian motion, continue on to Stochastic Processes 4: The Law of Large Numbers and the Central Limit Theorem.

The complete chain of statistical inference can be summarized as follows: the sampling design determines whether the data is representative, a statistic compresses the data into an interpretable number, an estimation method connects the statistic to the unknown parameter, bias and MSE assess the quality of a finite sample, and the law of large numbers together with the central limit theorem explain how an estimator stabilizes and fluctuates as the sample grows. Every step depends on explicit model assumptions; no formula on its own can compensate for a flawed sampling method or an inappropriate probability model.