"""Estimator quality metrics, weighted least squares, and the Cramer-Rao bound. The functions here are the working tools behind the estimation & detection arc: a way to score an estimator empirically (bias, variance, MSE), the weighted-least-squares fit that exploits a known noise model, and the Cramer-Rao lower bound that says how well any unbiased estimator could possibly do. For the linear Gaussian model the three meet: WLS attains the bound, so it is simultaneously the least-squares fit, the maximum likelihood estimate, and the minimum-variance unbiased estimator. """ import numpy as np def bias_variance_mse(estimates: np.ndarray, true_value: float | np.ndarray) -> tuple[float, float, float]: """Score an estimator from a Monte Carlo ensemble of its estimates. Parameters ---------- estimates : ndarray of shape (n_trials,) One estimate per independent realisation of the data. true_value : float The parameter value that generated the data. Returns ------- bias : float E[theta_hat] - theta, estimated by the sample mean. variance : float var(theta_hat), estimated with the 1/n (biased) normalisation so that the identity ``mse = bias**2 + variance`` holds exactly on the sample, not just in expectation. mse : float Mean squared error E[(theta_hat - theta)^2]. Notes ----- The decomposition MSE = bias^2 + variance is the reason "unbiased" is not the same as "good": an estimator may trade a little bias for a large reduction in variance and win on MSE. See the topic page for a worked case. """ estimates = np.asarray(estimates, dtype=float) if estimates.ndim != 1: raise ValueError("estimates must be a 1-D array of scalar estimates") if estimates.size < 2: raise ValueError("need at least two trials to estimate a variance") bias = float(estimates.mean() - true_value) variance = float(estimates.var()) # 1/n normalisation, see above mse = float(np.mean((estimates - true_value) ** 2)) return bias, variance, mse def weighted_least_squares(X: np.ndarray, y: np.ndarray, noise_var: np.ndarray | float | None = None ) -> tuple[np.ndarray, np.ndarray]: """Fit ``y = X @ theta + n`` with per-sample noise variances. Each observation is weighted by the inverse of its own noise variance, so precise samples pull harder on the fit than noisy ones. With ``noise_var=None`` (or a scalar) every sample is weighted equally and this reduces to ordinary least squares. Parameters ---------- X : ndarray of shape (n_samples, n_params) Design matrix: column j holds basis function j sampled at every observation point. y : ndarray of shape (n_samples,) Observations. noise_var : ndarray of shape (n_samples,), float, or None Known noise variance per sample. Must be strictly positive. Returns ------- theta : ndarray of shape (n_params,) The weighted least-squares estimate ``theta = (X^T W X)^-1 X^T W y`` with ``W = diag(1/noise_var)``. cov : ndarray of shape (n_params, n_params) Covariance of the estimate, ``(X^T W X)^-1``. When *noise_var* gives the true variances this is exact, and equals the Cramer-Rao bound for the Gaussian case: the estimator is efficient. With a scalar or None it is exact up to the unknown common noise variance factor. Raises ------ ValueError If the shapes disagree, a variance is non-positive, or the problem is underdetermined. """ X = np.atleast_2d(np.asarray(X, dtype=float)) y = np.asarray(y, dtype=float).ravel() n_samples, n_params = X.shape if y.size != n_samples: raise ValueError(f"X has {n_samples} rows but y has {y.size} entries") if n_samples < n_params: raise ValueError("underdetermined: fewer observations than parameters") if noise_var is None: w = np.ones(n_samples) else: var = np.broadcast_to(np.asarray(noise_var, dtype=float), (n_samples,)) if np.any(var <= 0): raise ValueError("noise_var must be strictly positive") w = 1.0 / var XtW = X.T * w # (n_params, n_samples) fisher = XtW @ X # X^T W X cov = np.linalg.inv(fisher) theta = cov @ (XtW @ y) return theta, cov def crlb_linear_gaussian(X: np.ndarray, noise_var: np.ndarray | float) -> np.ndarray: """Cramer-Rao lower bound for the linear Gaussian model. For ``y = X @ theta + n`` with ``n ~ N(0, diag(noise_var))`` and known noise variances, the Fisher information matrix is ``I = X^T W X`` with ``W = diag(1/noise_var)``, so no unbiased estimator of theta can have a covariance smaller than ``I^-1`` (in the positive-semidefinite sense). In particular the variance of any unbiased estimate of ``theta[j]`` is at least ``(I^-1)[j, j]``. Parameters ---------- X : ndarray of shape (n_samples, n_params) Design matrix, as in :func:`weighted_least_squares`. noise_var : ndarray of shape (n_samples,) or float Known noise variance per sample. Returns ------- crlb : ndarray of shape (n_params, n_params) The bound ``(X^T W X)^-1``. Its diagonal holds the per-parameter variance bounds. Notes ----- This is deliberately the same expression that :func:`weighted_least_squares` returns as its covariance. That is not a coincidence and not a shortcut: for this model WLS *attains* the bound, which is what makes it the minimum-variance unbiased estimator. The two functions are kept separate because they answer different questions: one fits data, the other bounds what any fit could achieve, and the second needs no data at all, only the measurement design. """ X = np.atleast_2d(np.asarray(X, dtype=float)) n_samples = X.shape[0] var = np.broadcast_to(np.asarray(noise_var, dtype=float), (n_samples,)) if np.any(var <= 0): raise ValueError("noise_var must be strictly positive") fisher = (X.T / var) @ X return np.linalg.inv(fisher) def polynomial_design(x: np.ndarray, powers) -> np.ndarray: """Design matrix whose columns are ``x**p`` for each p in *powers*. Convenience for fitting a polynomial with only *some* terms present, which is the point of the worked example on the topic page: knowing that the constant and quadratic terms are absent is prior knowledge, and prior knowledge is what buys accuracy. Parameters ---------- x : ndarray of shape (n_samples,) Sample points. powers : sequence of int Exponents to include, e.g. ``(1, 3, 4)`` for a polynomial with only the linear, cubic, and quartic terms. Returns ------- X : ndarray of shape (n_samples, len(powers)) """ x = np.asarray(x, dtype=float).ravel() powers = list(powers) if not powers: raise ValueError("need at least one power") return np.column_stack([x ** p for p in powers])