Statistics & Inference

Setting Up the Problem

Imagine you flip a coin 10 times and see 7 heads. What is the coin's bias? You cannot know for certain — you only have a sample. Estimation is the problem of inferring unknown quantities (called parameters) from data.

Formally, we assume data $x_1, x_2, \ldots, x_n$ are drawn i.i.d. (independently and identically distributed) from a probability distribution $P(x \mid \theta)$ parametrized by some unknown $\theta \in \Theta$. Our goal is to produce an estimator $\hat{\theta}(x_1, \ldots, x_n)$ — a function of the data — that is close to the true $\theta$.

The Likelihood Function

The likelihood of parameter $\theta$ given observations is: $$L(\theta) = P(x_1, \ldots, x_n \mid \theta) = \prod_{i=1}^n P(x_i \mid \theta)$$

Crucial distinction: The likelihood is NOT a probability distribution over $\theta$. It is a function of $\theta$ with the data fixed. We ask: for which $\theta$ was this particular dataset most probable to generate?

It is almost always easier to work with the log-likelihood $\ell(\theta) = \log L(\theta) = \sum_{i=1}^n \log P(x_i \mid \theta)$, which turns products into sums (and since $\log$ is monotone, maximizing $\ell$ is equivalent to maximizing $L$).

Evaluating Estimators

Two key properties:

• Bias: $\text{Bias}(\hat{\theta}) = E[\hat{\theta}] - \theta$. An estimator is unbiased if its expectation equals the true value.
• Variance: How much does $\hat{\theta}$ fluctuate across different datasets?

The bias-variance tradeoff is fundamental: $\text{MSE}(\hat{\theta}) = \text{Bias}^2 + \text{Variance}$. Sometimes a biased estimator with low variance beats an unbiased one with high variance.

import numpy as np

# Estimating the bias p of a coin from n flips
def sample_variance_unbiased(data):
    """Bessel's correction: divide by n-1, not n."""
    n = len(data)
    mean = np.mean(data)
    return np.sum((data - mean)**2) / (n - 1)  # unbiased estimator of sigma^2

def sample_variance_biased(data):
    """MLE for Gaussian variance — biased downward."""
    n = len(data)
    mean = np.mean(data)
    return np.sum((data - mean)**2) / n  # biased: E[s^2] = (n-1)/n * sigma^2

np.random.seed(42)
data = np.random.normal(loc=5.0, scale=2.0, n := 20, size=n)
print(f"True variance: 4.0")
print(f"Biased MLE:    {sample_variance_biased(data):.4f}")
print(f"Unbiased:      {sample_variance_unbiased(data):.4f}")

The MLE Principle

Maximum Likelihood Estimation finds the parameter value that maximizes the probability of observing the actual data: $$\hat{\theta}{\text{MLE}} = \arg\max{\theta} \ell(\theta) = \arg\max_{\theta} \sum_{i=1}^n \log P(x_i \mid \theta)$$

This is an intuitively compelling principle: choose the world that best explains what you saw.

MLE for Bernoulli (Coin Flips)

Let $x_i \in {0,1}$, $P(x_i=1 \mid p) = p$. Then: $$\ell(p) = \sum_i x_i \log p + (1-x_i)\log(1-p)$$

Setting $d\ell/dp = 0$: $$\hat{p}_{\text{MLE}} = \frac{#\text{heads}}{n}$$

The sample proportion is the MLE for Bernoulli probability — your intuition was right all along. MLE formalizes it.

MLE for Gaussian

For $x_i \sim \mathcal{N}(\mu, \sigma^2)$: $$\ell(\mu, \sigma^2) = -\frac{n}{2}\log(2\pi\sigma^2) - \frac{1}{2\sigma^2}\sum_i (x_i - \mu)^2$$

MLE gives $\hat{\mu} = \bar{x}$ (sample mean) and $\hat{\sigma}^2 = \frac{1}{n}\sum(x_i-\bar{x})^2$ — the biased sample variance. The unbiased correction ($n-1$) comes from a different principle (UMVUE theory).

Why MLE is So Powerful

MLE has remarkable asymptotic properties (as $n \to \infty$):

1. Consistency: $\hat{\theta}_{\text{MLE}} \to \theta$ in probability
2. Asymptotic normality: $\sqrt{n}(\hat{\theta}_{\text{MLE}} - \theta) \to \mathcal{N}(0, I(\theta)^{-1})$, where $I(\theta)$ is the Fisher Information
3. Efficiency: MLE achieves the Cramér-Rao lower bound asymptotically — no consistent estimator can have lower variance

The Fisher Information $I(\theta) = E\left[\left(\frac{d}{d\theta}\log P(x\mid\theta)\right)^2\right]$ measures how much information the data carries about $\theta$. High Fisher information = sharp likelihood peak = precise estimation.

import numpy as np
from scipy.optimize import minimize_scalar

# MLE for a Poisson distribution: P(x | lambda) = e^{-lambda} * lambda^x / x!
# Log-likelihood: sum_i [ x_i * log(lambda) - lambda - log(x_i!) ]
# MLE analytically: lambda_hat = sample_mean

def neg_log_likelihood_poisson(lam, data):
    if lam <= 0:
        return np.inf
    return -(np.sum(data * np.log(lam)) - len(data) * lam)

np.random.seed(0)
true_lambda = 3.5
data = np.random.poisson(true_lambda, size=100)

# Numerical MLE
result = minimize_scalar(neg_log_likelihood_poisson, bounds=(0.01, 20),
                         method='bounded', args=(data,))
print(f"True lambda:  {true_lambda}")
print(f"MLE (numeric):{result.x:.4f}")
print(f"Sample mean:  {data.mean():.4f}")  # Closed-form MLE

A Different Philosophy

In the frequentist view, $\theta$ is a fixed (unknown) constant and probability refers to long-run frequencies. MLE asks: what fixed $\theta$ maximizes the data probability?

In the Bayesian view, $\theta$ is itself a random variable with a prior distribution $P(\theta)$ encoding your beliefs before seeing data. After observing data $\mathbf{x}$, you update to the posterior:

$$P(\theta \mid \mathbf{x}) = \frac{P(\mathbf{x} \mid \theta) \cdot P(\theta)}{P(\mathbf{x})} \propto \underbrace{P(\mathbf{x} \mid \theta)}{\text{likelihood}} \cdot \underbrace{P(\theta)}{\text{prior}}$$

This is Bayes' theorem applied to inference. The denominator $P(\mathbf{x}) = \int P(\mathbf{x}\mid\theta)P(\theta)d\theta$ is the marginal likelihood (evidence) — a normalizing constant.

Conjugate Priors: Analytic Updates

For many likelihood/prior pairs, the posterior has the same family as the prior — these are conjugate priors.

Example: Beta-Binomial conjugacy.

• Prior: $p \sim \text{Beta}(\alpha, \beta)$
• Likelihood: $k$ heads in $n$ flips — Binomial$(n,p)$
• Posterior: $p \mid k \sim \text{Beta}(\alpha + k,\ \beta + n - k)$

Intuition: $\alpha$ and $\beta$ are like "pseudo-counts" of prior heads and tails. Evidence just adds real counts.

MAP vs Posterior Mean

• MAP (Maximum A Posteriori): $\hat{\theta}_{\text{MAP}} = \arg\max P(\theta\mid\mathbf{x})$ — the mode of the posterior. With uniform prior, MAP = MLE.
• Posterior Mean: $E[\theta \mid \mathbf{x}]$ — minimizes posterior expected squared error.

For Beta$(\alpha+k, \beta+n-k)$, the posterior mean is $\frac{\alpha+k}{\alpha+\beta+n}$ — a weighted blend of the prior belief $\alpha/(\alpha+\beta)$ and the data proportion $k/n$. As $n\to\infty$, the data drowns the prior.

Credible Intervals vs Confidence Intervals

A Bayesian 95% credible interval $[a,b]$ satisfies $P(a \le \theta \le b \mid \mathbf{x}) = 0.95$. This is what most people think a confidence interval means — but a confidence interval is a frequentist object with a subtly different interpretation (see below).

import numpy as np
from scipy.stats import beta
import matplotlib.pyplot as plt

# Bayesian coin-flip inference with Beta-Binomial conjugacy

alpha_prior, beta_prior = 2, 2  # weak prior: slightly favour fairness
np.random.seed(7)
true_p = 0.7

flips = np.random.binomial(1, true_p, size=30)
heads = flips.sum()
tails = len(flips) - heads

alpha_post = alpha_prior + heads
beta_post  = beta_prior  + tails

print(f"Observed: {heads} heads, {tails} tails")
print(f"Prior:     Beta({alpha_prior}, {beta_prior}) — mean={alpha_prior/(alpha_prior+beta_prior):.2f}")
print(f"Posterior: Beta({alpha_post}, {beta_post}) — mean={alpha_post/(alpha_post+beta_post):.2f}")

# 95% credible interval (highest-density interval approximated by quantiles)
lo, hi = beta.ppf([0.025, 0.975], alpha_post, beta_post)
print(f"95% Credible Interval: [{lo:.3f}, {hi:.3f}]")
print(f"True p = {true_p}")

The Most Important Theorem in Applied Statistics

Let $X_1, X_2, \ldots, X_n$ be i.i.d. random variables with mean $\mu$ and variance $\sigma^2 < \infty$. Define the sample mean $\bar{X}n = \frac{1}{n}\sum{i=1}^n X_i$. Then:

$$\sqrt{n}\left(\bar{X}_n - \mu\right) \xrightarrow{d} \mathcal{N}(0, \sigma^2) \quad \text{as } n \to \infty$$

Or equivalently, $\bar{X}_n \approx \mathcal{N}\left(\mu, \frac{\sigma^2}{n}\right)$ for large $n$.

This holds regardless of the underlying distribution of $X_i$ (provided $\sigma^2 < \infty$). Exponential, uniform, Poisson, Bernoulli — average them, and the distribution of the average becomes normal. The standard deviation of this distribution, $\sigma/\sqrt{n}$, is called the standard error of the mean.

Why Does This Happen?

The CLT emerges from the structure of the characteristic function (Fourier transform of the density). When you sum many independent small contributions, the characteristic function of the sum converges to that of a Gaussian: $e^{-t^2/2}$. The Gaussian is the attractor in the space of distributions under convolution-and-rescaling. This deep fact has analogues in physics (diffusion), information theory, and even random matrix theory.

When the CLT Fails

• Heavy tails: $\alpha$-stable distributions (e.g., Cauchy, $\sigma^2 = \infty$). The mean may not even exist. Averages converge to a stable law, not Gaussian.
• Dependence: The CLT generalizes to dependent sequences (martingale CLT, mixing conditions), but requires care.
• Very small $n$: The approximation needs roughly $n \ge 30$ for symmetric distributions, much more for skewed ones.

Practical Consequence

Because of the CLT, whenever we work with sample means (or sums), we can apply Gaussian statistics. This is why the $z$-test and $t$-test work even when individual measurements are not Gaussian.

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm

np.random.seed(42)
n_samples = 10000

# Try various underlying distributions
for dist_name, gen in [
    ("Uniform(0,1)",   lambda n: np.random.uniform(0, 1, n)),
    ("Exponential(1)", lambda n: np.random.exponential(1, n)),
    ("Bernoulli(0.3)", lambda n: np.random.binomial(1, 0.3, n)),
]:
    for n in [1, 5, 30, 100]:
        # Draw n_samples replications, each of size n
        means = np.array([gen(n).mean() for _ in range(n_samples)])
        print(f"{dist_name}, n={n:3d}: "
              f"mean={means.mean():.3f}, "
              f"std={means.std():.4f} "
              f"(theory: {1/np.sqrt(n)*0.5:.4f} approx)")

What a Confidence Interval Actually Means

A 95% confidence interval (CI) is a random interval — a function of the data — constructed such that if you repeated the experiment many times, 95% of the resulting intervals would contain the true parameter $\theta$.

Critical point: After you compute a specific interval $[2.3, 4.7]$, it is NOT correct to say "there is a 95% probability that $\theta \in [2.3, 4.7]$." The true $\theta$ is fixed; the interval is what's random. Either it contains $\theta$ (probability 1) or it doesn't (probability 0). The 95% refers to the procedure, not any particular realized interval.

Constructing a CI for the Mean

By the CLT, $\bar{X}_n \approx \mathcal{N}(\mu, \sigma^2/n)$. So: $$\frac{\bar{X}_n - \mu}{\sigma/\sqrt{n}} \approx \mathcal{N}(0,1)$$

Rearranging: a 95% CI is $$\bar{X}n \pm z{0.025} \cdot \frac{\sigma}{\sqrt{n}} = \bar{X}_n \pm 1.96 \cdot \frac{\sigma}{\sqrt{n}}$$

When $\sigma$ is unknown (almost always), we plug in the sample standard deviation $s$ and replace the $z$-quantile with a $t$-quantile from Student's $t$-distribution with $n-1$ degrees of freedom. For large $n$ the $t$ and $z$ quantiles are nearly identical.

The $t$-Distribution

Student's $t$ with $\nu$ degrees of freedom arises naturally from: $$T = \frac{\bar{X} - \mu}{s / \sqrt{n}} \sim t_{n-1}$$

It has heavier tails than the normal (more probability mass in extreme events), which accounts for our uncertainty about $\sigma$. As $\nu \to \infty$, $t_\nu \to \mathcal{N}(0,1)$.

Interpretation Analogy

Think of a CI construction as a fishing net. The net has a fixed mesh size (confidence level). You cast it many times; 95% of casts will land on the fish (the true parameter). After a single cast, the fish is either in your net or not — but you don't know which.

import numpy as np
from scipy import stats

def confidence_interval_mean(data, confidence=0.95):
    n = len(data)
    mean = np.mean(data)
    se = stats.sem(data)  # standard error = s / sqrt(n)
    # t-interval (works for any n; approaches z-interval as n -> inf)
    t_crit = stats.t.ppf((1 + confidence) / 2, df=n - 1)
    margin = t_crit * se
    return mean - margin, mean + margin

np.random.seed(0)
true_mean = 10.0

# Simulate coverage: how often does 95% CI contain true mean?
coverage_count = 0
for _ in range(10000):
    sample = np.random.exponential(scale=true_mean, size=30)  # NOT normal!
    lo, hi = confidence_interval_mean(sample, confidence=0.95)
    if lo <= true_mean <= hi:
        coverage_count += 1

print(f"Empirical coverage: {coverage_count/100:.1f}%")  # Should be ~95%

Framing a Scientific Question as a Binary Decision

Hypothesis testing is a formal decision procedure: given data, should we reject or retain a specific claim about the world?

• Null hypothesis $H_0$: the default, conservative claim (e.g., "the drug has no effect", "the coin is fair", "the two groups have equal means")
• Alternative hypothesis $H_1$: the claim we're trying to establish evidence for

The p-value

The p-value is the probability of observing data at least as extreme as what we saw, assuming $H_0$ is true: $$p = P(T \ge t_{\text{obs}} \mid H_0)$$

where $T$ is a test statistic — a summary of the data calibrated so that its distribution under $H_0$ is known.

The p-value is NOT: the probability that $H_0$ is true. That would require a prior (Bayesian approach). The p-value tells you how surprising the data are if the null were true.

Type I and Type II Errors

The significance level $\alpha$ (typically 0.05) is the maximum acceptable Type I error rate. Setting $\alpha = 0.05$ means: "I am willing to wrongly reject a true null 5% of the time."

The z-test and t-test

For testing $H_0: \mu = \mu_0$: $$z = \frac{\bar{X} - \mu_0}{\sigma/\sqrt{n}} \quad \text{(z-test, known } \sigma \text{)}$$ $$t = \frac{\bar{X} - \mu_0}{s/\sqrt{n}} \sim t_{n-1} \quad \text{(t-test, unknown } \sigma \text{)}$$

Multiple Testing

If you run 100 independent tests at $\alpha = 0.05$, you expect ~5 false positives by chance. The Bonferroni correction adjusts $\alpha' = \alpha/m$ for $m$ tests. The Benjamini-Hochberg procedure controls the False Discovery Rate (FDR) at level $q$, which is less conservative than Bonferroni and standard in genomics and ML.

import numpy as np
from scipy import stats

# One-sample t-test: is this population mean different from 0?
np.random.seed(1)
data = np.random.normal(loc=0.5, scale=2.0, size=40)  # true mean = 0.5

t_stat, p_value = stats.ttest_1samp(data, popmean=0)
print(f"t-statistic: {t_stat:.4f}")
print(f"p-value:     {p_value:.4f}")
print(f"Reject H0 at alpha=0.05: {p_value < 0.05}")

# Two-sample Welch t-test: do two groups differ?
group_a = np.random.normal(5.0, 1.5, size=50)
group_b = np.random.normal(5.8, 2.0, size=50)  # slightly higher mean
t2, p2 = stats.ttest_ind(group_a, group_b, equal_var=False)  # Welch
print(f"\nTwo-sample: t={t2:.3f}, p={p2:.4f}")

# Effect size: Cohen's d (what p-value CANNOT tell you)
d = (group_b.mean() - group_a.mean()) / np.sqrt(
    (group_a.std()**2 + group_b.std()**2) / 2
)
print(f"Cohen's d (effect size): {d:.3f}")

What Data Really Carries

A sufficient statistic $T(\mathbf{x})$ captures all the information in the data about parameter $\theta$. Formally, the conditional distribution $P(\mathbf{x} \mid T(\mathbf{x}), \theta)$ does not depend on $\theta$. Once you know $T$, the raw data carry no additional information about $\theta$.

The Fisher-Neyman factorization theorem characterizes sufficiency: $T$ is sufficient iff the likelihood factors as: $$P(\mathbf{x} \mid \theta) = h(\mathbf{x}) \cdot g(T(\mathbf{x}), \theta)$$

Example: For $X_i \sim \mathcal{N}(\mu, \sigma^2)$, the sufficient statistic for $\mu$ is $\bar{X}$ (and for $\sigma^2$, it is $\sum(X_i - \bar{X})^2$). This justifies compressing $n$ data points into 2 numbers.

The Exponential Family

An astonishing fraction of distributions used in practice belong to the exponential family: $$P(x \mid \boldsymbol{\eta}) = h(x) \exp\left(\boldsymbol{\eta}^\top T(x) - A(\boldsymbol{\eta})\right)$$

where $\boldsymbol{\eta}$ are natural parameters, $T(x)$ is the sufficient statistic, and $A(\boldsymbol{\eta}) = \log \int h(x)\exp(\boldsymbol{\eta}^\top T(x))dx$ is the log-partition function (ensures normalization).

Members: Gaussian, Bernoulli, Poisson, Exponential, Beta, Dirichlet, Gamma, Multinomial — essentially every named univariate distribution.

Why it matters for ML: The MLE for exponential family models always has the clean form "expected sufficient statistics equal observed sufficient statistics". Generalized Linear Models (logistic regression, Poisson regression) are all exponential family models with a linear predictor.

# The exponential family in action: logistic regression
# Bernoulli log-likelihood under sigmoid link

import numpy as np

def sigmoid(z):
    return 1 / (1 + np.exp(-z))

def log_likelihood_logistic(X, y, w):
    """Binary cross-entropy = negative Bernoulli log-likelihood."""
    probs = sigmoid(X @ w)
    # Clamp to avoid log(0)
    probs = np.clip(probs, 1e-12, 1 - 1e-12)
    return np.sum(y * np.log(probs) + (1 - y) * np.log(1 - probs))

def grad_log_likelihood(X, y, w):
    """Gradient: X^T (y - sigma(Xw)) — elegant because of exp family."""
    return X.T @ (y - sigmoid(X @ w))

# The gradient has the form: (observed sufficient stat) - (expected sufficient stat)
# This is the universal MLE condition for exponential families.

Bias-Variance in Estimation

Any estimator's mean squared error decomposes as: $$\text{MSE}(\hat{\theta}) = \underbrace{\left(E[\hat{\theta}] - \theta\right)^2}{\text{bias}^2} + \underbrace{\text{Var}(\hat{\theta})}{\text{variance}}$$

This is the estimation analogue of the ML generalization tradeoff. A more complex model (or estimator) can reduce bias but inflates variance, especially in small samples.

Example: Using the MLE $\hat{\sigma}^2 = \frac{1}{n}\sum(x_i-\bar{x})^2$ instead of the unbiased $s^2 = \frac{1}{n-1}\sum(x_i-\bar{x})^2$ introduces a small downward bias, but for $n \ge 30$ the MSE is actually lower for the biased MLE (less variance).

The Bootstrap

The bootstrap is a general-purpose technique for quantifying uncertainty in any estimator without relying on analytic formulas.

Algorithm (non-parametric bootstrap):

1. Have data $\mathbf{x} = {x_1, \ldots, x_n}$
2. For $b = 1, \ldots, B$:
   • Draw a bootstrap sample $\mathbf{x}^*$ of size $n$ with replacement from $\mathbf{x}$
   • Compute $\hat{\theta}^*(b)$ on the bootstrap sample
3. The bootstrap distribution of $\hat{\theta}^*$ approximates the sampling distribution of $\hat{\theta}$
4. Use percentiles of ${\hat{\theta}^*(b)}$ to form confidence intervals

Key insight: The relationship between $\hat{\theta}$ and the true $\theta$ is mimicked by the relationship between $\hat{\theta}^*$ and $\hat{\theta}$, because the empirical distribution $\hat{F}_n$ (which assigns $1/n$ to each data point) approximates the true $F$.

The bootstrap is invaluable when the estimator is complex (e.g., a median, a correlation coefficient, a neural network weight) and analytic variance formulas don't exist.

import numpy as np

def bootstrap_ci(data, statistic_fn, B=10000, alpha=0.05):
    """
    Non-parametric bootstrap confidence interval.
    statistic_fn: data -> scalar
    """
    n = len(data)
    boot_stats = np.empty(B)
    for b in range(B):
        sample = data[np.random.randint(0, n, size=n)]  # resample w/ replacement
        boot_stats[b] = statistic_fn(sample)
    lo = np.percentile(boot_stats, 100 * alpha / 2)
    hi = np.percentile(boot_stats, 100 * (1 - alpha / 2))
    return lo, hi, boot_stats

np.random.seed(99)
data = np.random.lognormal(mean=1.0, sigma=0.5, size=50)

# CI for the median (no simple analytic formula!)
lo, hi, _ = bootstrap_ci(data, np.median, B=5000)
print(f"Median estimate:  {np.median(data):.3f}")
print(f"95% Bootstrap CI: [{lo:.3f}, {hi:.3f}]")
print(f"True median:      {np.exp(1.0):.3f}")

A Floor on Estimation Uncertainty

Not all unbiased estimators are created equal. The Cramér-Rao Lower Bound (CRLB) gives the minimum variance that any unbiased estimator can achieve: $$\text{Var}(\hat{\theta}) \ge \frac{1}{n \cdot I(\theta)}$$

where the Fisher Information is: $$I(\theta) = E\left[\left(\frac{\partial}{\partial\theta}\log P(x\mid\theta)\right)^2\right] = -E\left[\frac{\partial^2}{\partial\theta^2}\log P(x\mid\theta)\right]$$

Intuition: The Fisher information measures the curvature of the log-likelihood. A sharply peaked likelihood (high curvature, high information) means the data strongly constrains $\theta$, so estimation can be precise. A flat likelihood means $\theta$ is hard to identify — lower bound on variance is large.

Example: Bernoulli

$P(x\mid p) = p^x(1-p)^{1-x}$, $\log P = x\log p + (1-x)\log(1-p)$.

$\frac{d}{dp}\log P = \frac{x}{p} - \frac{1-x}{1-p}$, so $I(p) = \frac{1}{p(1-p)}$.

CRLB: $\text{Var}(\hat{p}) \ge \frac{p(1-p)}{n}$.

The sample mean achieves this: $\text{Var}(\bar{X}) = p(1-p)/n$. It is efficient — no unbiased estimator can do better.

Consequences

• An estimator achieving the CRLB is called efficient (or UMVUE — Uniformly Minimum Variance Unbiased Estimator)
• MLE is asymptotically efficient (achieves CRLB as $n\to\infty$)
• The CRLB is fundamental in signal processing, where it bounds how precisely you can estimate a signal parameter from noisy measurements

# Verifying the CRLB numerically for Bernoulli MLE
import numpy as np

np.random.seed(42)
true_p = 0.4
n = 100
M = 50000  # number of Monte Carlo repetitions

# Simulate M datasets of size n; compute MLE (= sample mean) for each
estimates = np.array([
    np.random.binomial(1, true_p, size=n).mean()
    for _ in range(M)
])

emp_var = np.var(estimates)
crlb    = true_p * (1 - true_p) / n

print(f"Empirical variance of MLE: {emp_var:.6f}")
print(f"Cramér-Rao lower bound:    {crlb:.6f}")
print(f"Ratio (should be ~1):      {emp_var/crlb:.4f}")
# MLE is efficient — it exactly achieves the CRLB