06 - Statistical Inference: Estimators and Properties

You have measured the same quantity twenty times and you want to say how much the measurements spread out. The sample variance formula has a denominator: the textbook says 19, while the calculator has two keys, one that divides by 20 and one by 19. The short answer is that neither is “the right one”: with 19 the estimator is unbiased, with 20 it coincides with the maximum likelihood one, and with 21 it has the smallest mean squared error of the three, as I compute below. Choosing means deciding which property you care about.

This article, the sixth in the Probability and Statistics series, gives you the vocabulary to do it: what an estimator is, how its quality is measured, how well one can do at best, and how a good estimator is built. You need random variables, expected value and variance, and it helps to have read the article on the central limit theorem and the law of large numbers.

Parameter, estimator, estimate

Imagine weighing the same object twenty times on an imprecise scale. Assume that each weighing X_i is a Normal random variable with mean μ (the true weight) and variance σ² (the imprecision of the scale), and that the weighings are independent. The numbers μ and σ² are the parameters: they describe the distribution and you do not know them. A statistic is any function of the twenty weighings, such as the mean or the largest one. An estimator is a statistic chosen to get close to a parameter; the estimate is the number the estimator takes on the sample you actually hold.

The distinction matters because the estimator, before the data are collected, is a random variable: repeating the twenty weighings you would get a different mean each time. Its distribution is called the sampling distribution, and it is this, not the single estimate, that lets you judge the estimator.

If I call T an estimator of a parameter θ, its quality is described by three numbers: the bias b = E[T] − θ, that is, by how much it is wrong on average; the variance Var(T), that is, how much it fluctuates from sample to sample; and the mean squared error (MSE) E[(T − θ)²], which measures the average squared distance from the target.

The three numbers are tied by an identity: MSE = Var(T) + b². The proof takes one line: write T − θ as (T − E[T]) + (E[T] − θ), square it and take the expected value; the cross product is zero because E[T − E[T]] = 0. The identity says that an estimator can trade a little bias for a smaller variance, and this is why “unbiased” does not mean “best”. The identity, under the name “quadratic risk = bias² + variance”, is in the parametric inference lecture of the MIT 18.650 course.

Dividing by n−1, by n or by n+1

For a Normal with n observations, let Q = Σ(X_i − X̄)² be the sum of the squared deviations from the sample mean. A classical result for Normal samples says that Q/σ² follows a chi-square distribution with n − 1 degrees of freedom, so it has mean n − 1 and variance 2(n − 1). Everything follows from here, for the estimator Q/d with divisor d: the mean is σ²(n − 1)/d and the variance is 2σ⁴(n − 1)/d².

With d = n − 1 the bias is zero; with d = n the mean is σ²(n − 1)/n, so the estimator always underestimates σ² by σ²/n. The intuitive reason is that the deviations are measured from the mean of the sample, which is the closest possible point to the data, and not from the true mean: Q is by construction smaller than it would be if measured from μ.

The program checks these formulas with a simulation: it generates 200,000 samples of n = 20 Normal weighings with σ² = 4, computes Q for each one and compares, for the three divisors, the mean and the MSE obtained from the simulation with the exact ones. The fixed seed makes the run repeatable, and the standard error of the mean of the squared errors is printed alongside, because a simulation only makes sense if you state how far it can be off.

# Three variance estimators: divisor n-1, n, n+1 (Normal sample, n = 20)
import numpy as np

n, sigma2 = 20, 4.0
R = 200_000                      # simulated samples
rng = np.random.default_rng(42)
X = rng.normal(0.0, np.sqrt(sigma2), size=(R, n))
somma_quadrati = ((X - X.mean(axis=1, keepdims=True)) ** 2).sum(axis=1)

print(f"n = {n}, sigma^2 = {sigma2}, simulated samples = {R}")
print("divisor  exact mean  sim. mean  exact MSE  sim. MSE (+/- std. error)")
for d in (n - 1, n, n + 1):
    stima = somma_quadrati / d
    # Exact: somma_quadrati / sigma^2 ~ chi-square with n-1 degrees of freedom,
    # so mean = sigma^2 (n-1)/d and variance = 2 sigma^4 (n-1)/d^2
    media = sigma2 * (n - 1) / d
    var = 2 * sigma2**2 * (n - 1) / d**2
    mse = var + (media - sigma2) ** 2
    err2 = (stima - sigma2) ** 2
    print(f"{d:7d}  {media:10.4f}  {stima.mean():9.4f}  {mse:9.4f}  "
          f"{err2.mean():8.4f} +/- {err2.std() / np.sqrt(R):.4f}")

# Output:
# n = 20, sigma^2 = 4.0, simulated samples = 200000
# divisor  exact mean  sim. mean  exact MSE  sim. MSE (+/- std. error)
#      19      4.0000     3.9999     1.6842    1.6896 +/- 0.0061
#      20      3.8000     3.7999     1.5600    1.5649 +/- 0.0052
#      21      3.6190     3.6189     1.5238    1.5283 +/- 0.0046

Read the output line by line. With divisor 19 the mean is 4.0000, that is, there is no bias, but the MSE is 1.6842. With divisor 20 the mean drops to 3.8, the bias is −0.2 and its square is 0.04, but the variance goes from 1.6842 to 1.52 and the overall MSE falls to 1.56. With 21 the bias is −0.381, the variance is 1.3787 and the MSE reaches 1.5238, the minimum of the three.

The simulated columns fall within about one standard error of the exact values. The gain is small (about 10% of MSE between 19 and 21) and shrinks as n grows; this is why in practice the divisor 19 is used because it is unbiased, and the variances of different groups can be averaged without accumulating a systematic error, not because it is the most precise.

What happens as n grows

A good estimator improves with more data. It is called consistent if, as n tends to infinity, it converges in probability to the parameter: the probability of being off by more than a fixed tolerance tends to zero. The sample mean is consistent by the law of large numbers. A useful counterexample is the estimator “take only the first weighing, X₁”: it is unbiased, but its variance stays σ² whatever n is, so it converges to nothing. Consistency and absence of bias are independent properties: the estimator with divisor n is biased and consistent, X₁ is unbiased and not consistent.

Consistency says where the estimator goes, not how fast it gets close. This is why asymptotic normality is also studied: for many estimators √n(T_n − θ) behaves, for large n, like a Normal with mean zero and variance V(θ).

For the maximum likelihood estimator, under regularity conditions, the limiting variance is the inverse of the Fisher information, which I discuss shortly: the theorem with its conditions (identifiable model, support that does not depend on θ, parameter not on the boundary of the parameter space) is stated in the maximum likelihood lecture of the MIT 18.650 course. Asymptotic normality is what makes it possible to build confidence intervals, as you will see in the article on hypothesis tests.

There are four properties asked of an estimator, and each answers a different question:

• Unbiased: on average, over many samples, does it hit the parameter?
• Consistent: with more data, does it get closer to the parameter?
• Efficient: does it have the smallest possible variance among those of its class?
• Sufficient: does it retain all the information about the parameter contained in the sample?

None of the four implies the others. The first two concern average behaviour and behaviour at the limit, the third compares estimators with each other, the fourth concerns information and not error. The rest of the article shows them one at a time on an example you can redo, starting from the lower bound on the variance, which is what gives a precise meaning to the word “efficient”.

Fisher information and the Cramér-Rao bound

How much information about θ does one observation carry? You answer by looking at the log-likelihood, that is, the logarithm of the density (or of the probability) seen as a function of θ. Its derivative with respect to θ is called the score, u(x, θ). The Fisher information is I(θ) = Var(u) = −E[u′]: if the log-likelihood is peaked (the derivative changes quickly), the data distinguish a value of θ well from its neighbours and I is large. The definition, with the equivalence between the two forms, is in the same MIT lecture on maximum likelihood.

For a Bernoulli with parameter p, the log-probability of x ∈ {0, 1} is x log p + (1 − x) log(1 − p), the score is x/p − (1 − x)/(1 − p), and −u′ = x/p² + (1 − x)/(1 − p)². The expected value is 1/p + 1/(1 − p) = 1/(p(1 − p)). The information is minimal at p = 1/2, where the coin is most unpredictable, and grows without limit towards 0 and 1.

The Cramér-Rao bound links information and variance: for every unbiased estimator of θ based on n independent observations, and under regularity conditions, Var(T) ≥ 1/(n I(θ)). The simplest condition to check is that the support of the distribution does not depend on θ: the Uniform on [0, θ] does not satisfy it, and for it the bound does not hold. For the Bernoulli the bound is p(1 − p)/n, which is exactly the variance of the sample frequency: that estimator is efficient.

For the exponential distribution with parameter λ (the rate, with mean 1/λ) the information is I(λ) = 1/λ², and the bound is λ²/n. Here the picture is less tidy than in the more hurried textbooks. I call T the sum of the n observations; T follows a Gamma distribution and 1/T an inverse Gamma, whose exact moments are returned by scipy.stats.invgamma (the SciPy documentation gives its definition).

Estimators of the form c/T include the maximum likelihood one (c = n, that is, 1/X̄) and the unbiased one (c = n − 1). The program computes for both, and for c = n − 2, the exact mean, variance and MSE, and compares them with a simulation.

# Cramer-Rao for Exp(lambda): the bound is not reached at finite n, not even by the unbiased estimator
import numpy as np
from scipy import stats

lam, n, R = 2.5, 20, 200_000
rng = np.random.default_rng(100)
X = rng.exponential(scale=1 / lam, size=(R, n))
T = X.sum(axis=1)                      # sufficient statistic: T ~ Gamma(n, scale 1/lam)

# Fisher information of one observation: variance of the score u(x) = 1/lam - x
I_sim = (1 / lam - X[:, 0]).var()
print(f"I(lambda): theoretical {1 / lam**2:.4f}, from the simulated score {I_sim:.4f}")
bound = 1 / (n * (1 / lam**2))
print(f"Cramer-Rao bound for unbiased estimators: {bound:.4f}")

# 1/T follows an inverse Gamma(n, scale lam): exact moments with scipy
inv = stats.invgamma(a=n, scale=lam)
m1, m2 = inv.mean(), inv.var() + inv.mean() ** 2
print("estimator c/T        exact mean  exact var  exact MSE  sim. var   sim. mean")
for c, nome in ((n, "MLE c=n"), (n - 1, "unbiased c=n-1"), (n - 2, "minimum MSE c=n-2")):
    media = c * m1
    var = c**2 * (m2 - m1**2)
    mse = var + (media - lam) ** 2
    s = c / T
    print(f"{nome:20s} {media:9.4f} {var:11.4f} {mse:11.4f} {s.var():9.4f} {s.mean():10.4f}")

# Output:
# I(lambda): theoretical 0.1600, from the simulated score 0.1586
# Cramer-Rao bound for unbiased estimators: 0.3125
# estimator c/T        exact mean  exact var  exact MSE  sim. var   sim. mean
# MLE c=n                 2.6316      0.3847      0.4020    0.3862     2.6304
# unbiased c=n-1          2.5000      0.3472      0.3472    0.3485     2.4988
# minimum MSE c=n-2       2.3684      0.3116      0.3289    0.3128     2.3673

The bound for n = 20 and λ = 2.5 is 0.3125. The maximum likelihood estimator has mean 2.6316, hence a bias of +0.1316, and variance 0.3847: the bound does not concern it, because it is not unbiased. The unbiased one (c = 19) has variance 0.3472, higher than the bound: it is the best among the unbiased ones, but at finite n it does not reach the bound; the variance is λ²/(n − 2), and the ratio to the bound, n/(n − 2), tends to 1 only for large n.

With c = 18 the variance (0.3116) falls even below the bound, and there is no contradiction, because the estimator is biased, and its MSE (0.3289) is the lowest of the three. The “I(lambda)” line also shows the simulation of the score: 0.1586 against a theoretical 0.16, with a standard error of about 0.001.

From information to the confidence interval

If the maximum likelihood estimator is approximately Normal with variance 1/(n I(θ)), then an interval going from T − z/√(n I) to T + z/√(n I), with z = 1.96, should contain θ in 95% of the samples. With I estimated from the data it is the Wald interval.

For a proportion p it becomes f ± 1.96 √(f(1 − f)/n), with f the sample frequency. “Should” is the key word: the program computes the true coverage, that is, the share of samples for which the interval actually contains p = 0.10, in two ways. With the exact formula, summing the binomial probabilities of the values of k for which the interval contains p, and with 200,000 simulated samples for each n.

# Actual coverage of the 95% Wald interval for a proportion p (Bernoulli)
import numpy as np
from scipy import stats

p, R, z = 0.10, 200_000, stats.norm.ppf(0.975)
rng = np.random.default_rng(2024)
print(f"p = {p}, z = {z:.4f}, samples per row = {R}")
print("   n  exact coverage  simulated coverage (+/- std. err.)")
for n in (20, 50, 100, 500, 2000):
    k = rng.binomial(n, p, size=R)          # number of successes in each sample
    ph = k / n
    se = np.sqrt(ph * (1 - ph) / n)         # 1/sqrt(n I(p)) with I estimated at ph
    copre = ((ph - z * se <= p) & (p <= ph + z * se)).mean()
    # Exact: sum of the binomial probabilities of the k for which the interval contains p
    kk = np.arange(n + 1)
    pk = kk / n
    sek = np.sqrt(pk * (1 - pk) / n)
    contiene = (pk - z * sek <= p) & (p <= pk + z * sek)
    esatta = stats.binom.pmf(kk, n, p)[contiene].sum()
    err = np.sqrt(copre * (1 - copre) / R)
    print(f"{n:4d}  {esatta:14.4f}  {copre:12.4f} +/- {err:.4f}")

# Output:
# p = 0.1, z = 1.9600, samples per row = 200000
#    n  exact coverage  simulated coverage (+/- std. err.)
#   20          0.8760        0.8766 +/- 0.0007
#   50          0.8789        0.8783 +/- 0.0007
#  100          0.9324        0.9325 +/- 0.0006
#  500          0.9427        0.9433 +/- 0.0005
# 2000          0.9462        0.9463 +/- 0.0005

With twenty observations the “95%” interval covers the parameter in 87.6% of the samples, and with fifty in 87.9%; about 2000 observations are needed to get close to 95% (94.6%). The simulation confirms the exact formula within a standard error of 0.0005-0.0007. An asymptotic property, that is, one true in the limit, guarantees nothing for a finite n: when p is close to 0 the Normal is a poor model of the distribution of f, and when k = 0 the interval collapses to a point that never contains p.

The correct reading of the “95%” is this: it is the frequency with which the procedure contains the parameter in the long run, and the table shows that sometimes it is much lower than stated. Once computed, the interval either contains p or it does not: it makes no sense to say that it contains it “with probability 95%”.

Sufficiency and Rao-Blackwell

A statistic T is sufficient for θ if, once the value of T is known, the distribution of the sample no longer depends on θ: the rest of the data is noise that says nothing about the parameter. For a Bernoulli sample, the sum of the successes T = ΣX_i is sufficient, because given the sum all the arrangements of the successes are equally likely whatever p is. It is recognised with the factorization theorem: T is sufficient when the joint density can be written as the product of a term that depends on the data only through T (and on θ) and one that does not depend on θ.

The Rao-Blackwell theorem gives sufficiency a practical use: if W is any unbiased estimator and T is sufficient, the estimator E[W | T] is still unbiased and has smaller or equal variance. Conditioning on T amounts to averaging over all the samples that have the same value of T, and averaging cannot increase the variance.

In the Bernoulli example, you start from W = X₁ (the first coin alone) and condition on the sum: by symmetry E[X₁ | T = t] = t/n, so the improved estimator is the sample frequency. The program checks it without trusting the symmetry: among the samples with sum 3, 6 and 10 it computes the mean of X₁ and compares it with t/n.

# Rao-Blackwell on Bernoulli(p): from X_1 to the sample mean, conditioning on the sum
import numpy as np

p, n, R = 0.3, 20, 1_000_000
rng = np.random.default_rng(3)
X = (rng.random(size=(R, n)) < p)
W = X[:, 0].astype(float)        # unbiased but noisy estimator: uses a single observation
T = X.sum(axis=1)                # sufficient statistic: number of successes
RB = T / n                       # E[W | T] = T/n

print(f"W = X_1     : mean {W.mean():.4f} (+/- {W.std() / np.sqrt(R):.4f}), variance {W.var():.4f} (exact {p * (1 - p):.4f})")
print(f"E[W|T] = T/n: mean {RB.mean():.4f} (+/- {RB.std() / np.sqrt(R):.4f}), variance {RB.var():.5f} (exact {p * (1 - p) / n:.5f})")
# Direct check of E[X_1 | T = t] = t/n, on the samples with sum t
for t in (3, 6, 10):
    sel = T == t
    err = W[sel].std() / np.sqrt(sel.sum())
    print(f"T = {t:2d}: samples {sel.sum():6d}, mean of X_1 {W[sel].mean():.4f} (+/- {err:.4f}), t/n = {t / n:.4f}")

# Output:
# W = X_1     : mean 0.3000 (+/- 0.0005), variance 0.2100 (exact 0.2100)
# E[W|T] = T/n: mean 0.2999 (+/- 0.0001), variance 0.01051 (exact 0.01050)
# T =  3: samples  71963, mean of X_1 0.1490 (+/- 0.0013), t/n = 0.1500
# T =  6: samples 191630, mean of X_1 0.2996 (+/- 0.0010), t/n = 0.3000
# T = 10: samples  30833, mean of X_1 0.5016 (+/- 0.0028), t/n = 0.5000

The variance of X₁ is 0.21 and that of T/n is 0.0105: exactly the factor n = 20, as predicted by the formulas p(1 − p) and p(1 − p)/n. Both have mean 0.3 within the printed error. The last three lines confirm the conditional identity: 0.1490 against 0.15, 0.2996 against 0.30, 0.5016 against 0.50, with standard errors of 0.001-0.003. Out of the twenty data points, the sum alone is enough to obtain the best estimator; the position of the successes adds nothing.

Worked exercise: the probability of zero events in a Poisson process

The number of failures per day of a machine is, by assumption, Poisson with parameter λ, and you want to estimate the probability of a day without failures, which is e^−λ. You have observed n days. A first unbiased estimator is W = 1 if there are no failures on the first day and 0 otherwise: it has expected value e^−λ, but it uses a single day.

The sum T of the failures over the n days is sufficient and follows a Poisson with parameter nλ. Given T = t, the t failures are distributed over the n days like t independent throws, each landing on day 1 with probability 1/n; the probability that day 1 stays empty is (1 − 1/n)^t. Rao-Blackwell therefore gives the estimator (1 − 1/n)^T.

Since the Poisson family is an exponential family, T is also complete (no non-trivial function of T has zero expected value for every λ), and the Lehmann-Scheffé theorem guarantees that this estimator is the unique unbiased one of minimum variance. Its variance is obtained from the generating function of a Poisson: with a = 1 − 1/n, E[a^2T] = exp(nλ(a² − 1)), hence Var = exp(nλ(a² − 1)) − e^−2λ. The program compares this estimator, its exact variance and the “natural” estimator e^−X̄ (λ is replaced with the sample mean), for λ = 2 and n = 10.

# Worked exercise: UMVUE of exp(-lambda) for a Poisson(lambda) sample, via Rao-Blackwell
import numpy as np

lam, n, R = 2.0, 10, 200_000
rng = np.random.default_rng(11)
X = rng.poisson(lam, size=(R, n))
T = X.sum(axis=1)                       # sufficient and complete: T ~ Poisson(n lam)

W = (X[:, 0] == 0).astype(float)        # unbiased estimator of P(X=0) = exp(-lam)
a = 1 - 1 / n
RB = a ** T                             # E[W | T = t] = (1 - 1/n)^t
plug = np.exp(-X.mean(axis=1))          # plug-in estimator exp(-mean), for comparison

target = np.exp(-lam)
var_w = target * (1 - target)
var_rb = np.exp(n * lam * (a**2 - 1)) - target**2     # from the Poisson(n lam) generating function
print(f"value to estimate exp(-lambda) = {target:.5f}")
print(f"W = 1{{X_1=0}}   : mean {W.mean():.5f}, variance {W.var():.5f} (exact {var_w:.5f})")
print(f"(1-1/n)^T      : mean {RB.mean():.5f}, variance {RB.var():.6f} (exact {var_rb:.6f})")
print(f"exp(-mean)     : mean {plug.mean():.5f}, variance {plug.var():.6f}, MSE {((plug - target) ** 2).mean():.6f}")

# Output:
# value to estimate exp(-lambda) = 0.13534
# W = 1{X_1=0}   : mean 0.13445, variance 0.11638 (exact 0.11702)
# (1-1/n)^T      : mean 0.13521, variance 0.004050 (exact 0.004055)
# exp(-mean)     : mean 0.14895, variance 0.004407, MSE 0.004593

The target is e⁻² = 0.13534. The raw estimator W has variance 0.117 (exact 0.11702), enormous compared with the value to be estimated; the Rao-Blackwell estimator has variance 0.004055 (simulated 0.004050), about 29 times smaller, and mean 0.13521, consistent with the unbiased one.

The plug-in estimator e^−X̄, the first one that would come to mind, is biased: its simulated mean is 0.14895, an overestimate of about 0.0136 due to the fact that e^−x is a convex function (Jensen’s inequality), and it has MSE 0.004593, worse than that of the unbiased estimator (0.004055). Here, unlike the case of the variance, the unique unbiased estimator of minimum variance wins on both fronts.

Moments versus maximum likelihood

The method of moments builds an estimator by equating the theoretical moments, which depend on θ, to the moments computed on the sample, and solving for θ; the justification of this procedure is in the chapter of the MIT course on the method of moments. For a Gamma with shape α and scale s, the mean and variance are αs and αs²; solving, α = mean²/variance, and it is enough to substitute the sample values.

Maximum likelihood instead chooses the parameter that makes the observed data most probable, that is, it maximizes the sum of the logarithms of the densities. For the Gamma there is no closed formula and a numerical optimizer is needed: scipy.stats.gamma.fit(x, floc=0) does it, and floc=0 fixes the location parameter at zero, as provided by the documentation of fit; without that argument the estimate would have one more parameter.

The program generates 2000 samples of n = 100 from a Gamma with α = 3 and scale 2, estimates α both ways, and compares the variances with the asymptotic Cramér-Rao bound, which is obtained by inverting the two-parameter Fisher information matrix, whose (α, α) entry requires the second derivative of the logarithm of the gamma function, the trigamma scipy.special.polygamma(1, α). The bound is an asymptotic reference, not a tight bound at n = 100.

# Gamma(alpha, scale): method of moments versus maximum likelihood for the shape alpha
import numpy as np
from scipy import stats
from scipy.special import polygamma

alfa, scala, n, R = 3.0, 2.0, 100, 2000
rng = np.random.default_rng(50)
mm = np.empty(R)
ml = np.empty(R)
for r in range(R):
    x = rng.gamma(alfa, scala, size=n)
    mm[r] = x.mean() ** 2 / x.var(ddof=1)            # moments: mean = alfa*scala, variance = alfa*scala^2
    ml[r] = stats.gamma.fit(x, floc=0)[0]            # maximum likelihood, location fixed at 0

# Asymptotic bound: (I^-1)_{alfa,alfa}/n with I = [[psi'(alfa), 1/scala], [1/scala, alfa/scala^2]]
I = np.array([[polygamma(1, alfa), 1 / scala], [1 / scala, alfa / scala**2]])
cr = np.linalg.inv(I)[0, 0] / n
print(f"true alpha = {alfa}, n = {n}, repetitions = {R}")
print(f"asymptotic Cramer-Rao bound for alpha: variance {cr:.4f}")
print(f"moments : mean {mm.mean():.4f}, variance {mm.var():.4f}, ratio to the bound {mm.var() / cr:.3f}")
print(f"MLE     : mean {ml.mean():.4f}, variance {ml.var():.4f}, ratio to the bound {ml.var() / cr:.3f}")

# Output:
# true alpha = 3.0, n = 100, repetitions = 2000
# asymptotic Cramer-Rao bound for alpha: variance 0.1623
# moments : mean 3.0932, variance 0.2603, ratio to the bound 1.604
# MLE     : mean 3.0841, variance 0.1905, ratio to the bound 1.173

The asymptotic bound for the variance of α is 0.1623. The method of moments has variance 0.2603, about 1.6 times the bound; maximum likelihood 0.1905, about 1.17 times. With 2000 repetitions a variance is known with a relative error of about 3%, so the difference between the two is well outside the noise.

Both have a small average overestimate (3.09 and 3.08 against 3). The price of the gain is an optimizer in place of a formula, with its possible non-convergence. In practice the two are used together: the moments give a starting point and a sanity check, the likelihood the final estimate.

Three mistakes that are really made

• “Unbiased means better”: it is false already for the variance of a Normal, where the divisor 21 has the lowest MSE.
• “No estimator can have variance below the Cramér-Rao bound”: the bound holds for unbiased estimators of a regular model, and in the exponential example a biased estimator exceeds it.
• “A 95% interval contains the parameter with probability 95%”: the 95% describes the procedure over many samples, not the single interval already computed, and at finite n it can be much less than 95%.

The three mistakes have the same origin: using an asymptotic or average property as if it were a guarantee on the single case. Each one is dismantled by a program of a few lines, and the simulation, when there is an exact formula to compare it with, is a cheap way to notice the mistake before writing it in a report.

You can now tell an estimator from an estimate, compute the bias, variance and MSE of an estimator (and choose between two by looking at the MSE), compute the Fisher information and the Cramér-Rao bound of a simple model, apply Rao-Blackwell to a raw estimator, and check with a simulation whether a confidence interval covers as much as it claims. In the next article you will use these tools to decide with a hypothesis test, and in the article on ANOVA you will extend them to more groups.

The same series also includes Discrete Random Variables, Multivariate Distributions and Covariance, Simple and Multiple Linear Regression, Bayesian Estimation: Prior, Posterior, and MAP, Bootstrap and Resampling Methods and Exploratory Data Analysis (EDA) with Python.