Static Estimation

Static parameter estimation algorithms.

Static estimation module.

This module provides methods for static (batch) estimation including least squares variants and robust estimation techniques.

class pytcl.static_estimation.LSResult(x, residuals, rank, singular_values)

Bases: NamedTuple

Result of least squares estimation.

x

Estimated parameters.

Type:

ndarray

residuals

Residual vector (y - A @ x).

Type:

ndarray

rank

Effective rank of the design matrix.

Type:

int

singular_values

Singular values of the design matrix.

Type:

ndarray

x: ndarray[tuple[Any, ...], dtype[floating]]

Alias for field number 0

residuals: ndarray[tuple[Any, ...], dtype[floating]]

Alias for field number 1

rank: int

Alias for field number 2

singular_values: ndarray[tuple[Any, ...], dtype[floating]]

Alias for field number 3

class pytcl.static_estimation.WLSResult(x, residuals, covariance, weighted_residual_sum)

Bases: NamedTuple

Result of weighted least squares estimation.

x

Estimated parameters.

Type:

ndarray

residuals

Residual vector.

Type:

ndarray

covariance

Estimated covariance of parameters.

Type:

ndarray

weighted_residual_sum

Sum of weighted squared residuals.

Type:

float

x: ndarray[tuple[Any, ...], dtype[floating]]

Alias for field number 0

residuals: ndarray[tuple[Any, ...], dtype[floating]]

Alias for field number 1

covariance: ndarray[tuple[Any, ...], dtype[floating]]

Alias for field number 2

weighted_residual_sum: float

Alias for field number 3

class pytcl.static_estimation.TLSResult(x, residuals_A, residuals_b, rank)

Bases: NamedTuple

Result of total least squares estimation.

x

Estimated parameters.

Type:

ndarray

residuals_A

Corrections to the design matrix A.

Type:

ndarray

residuals_b

Corrections to the observation vector b.

Type:

ndarray

rank

Effective rank used in the solution.

Type:

int

x: ndarray[tuple[Any, ...], dtype[floating]]

Alias for field number 0

residuals_A: ndarray[tuple[Any, ...], dtype[floating]]

Alias for field number 1

residuals_b: ndarray[tuple[Any, ...], dtype[floating]]

Alias for field number 2

rank: int

Alias for field number 3

pytcl.static_estimation.ordinary_least_squares(A, b, rcond=None)

Ordinary least squares estimation.

Solves the linear least squares problem:

min_x ||A @ x - b||_2

Parameters:

• A (array_like) – Design matrix of shape (m, n) where m >= n.

• b (array_like) – Observation vector of shape (m,) or (m, k).

• rcond (float, optional) – Cutoff for small singular values. Values smaller than rcond * largest_singular_value are treated as zero. Default is machine precision * max(m, n).

Returns:

result – Named tuple containing: - x: Estimated parameters of shape (n,) or (n, k) - residuals: Residual vector - rank: Effective rank of A - singular_values: Singular values of A

Return type:

LSResult

Examples

>>> A = np.array([[1, 1], [1, 2], [1, 3]])
>>> b = np.array([1, 2, 2])
>>> result = ordinary_least_squares(A, b)
>>> result.x  # Fitted line parameters
array([0.66666667, 0.5       ])

Notes

Uses SVD-based solution for numerical stability.

pytcl.static_estimation.weighted_least_squares(A, b, W=None, weights=None)

Weighted least squares estimation.

Solves the weighted linear least squares problem:

min_x (b - A @ x)^T W (b - A @ x)

Parameters:

• A (array_like) – Design matrix of shape (m, n).

• b (array_like) – Observation vector of shape (m,).

• W (array_like, optional) – Weight matrix of shape (m, m). If provided, weights is ignored. Should be positive definite.

• weights (array_like, optional) – Diagonal weights of shape (m,). Used only if W is None. Equivalent to W = diag(weights).

Returns:

result – Named tuple containing: - x: Estimated parameters of shape (n,) - residuals: Residual vector - covariance: Estimated parameter covariance - weighted_residual_sum: Weighted sum of squared residuals

Return type:

WLSResult

Examples

>>> A = np.array([[1, 1], [1, 2], [1, 3]])
>>> b = np.array([1, 2, 2])
>>> weights = np.array([1, 2, 1])  # Higher weight on middle observation
>>> result = weighted_least_squares(A, b, weights=weights)

Notes

The covariance of the estimated parameters is (A^T W A)^{-1}.

For measurement noise with covariance R, use W = R^{-1} to get the minimum variance unbiased estimator (MVUE).

pytcl.static_estimation.total_least_squares(A, b, rank=None)

Total least squares (TLS) estimation.

Solves the errors-in-variables problem where both the design matrix and observations may have errors:

min_{E, r} ||[E | r]||_F subject to (A + E) @ x = b + r

Parameters:

• A (array_like) – Design matrix of shape (m, n) where m >= n.

• b (array_like) – Observation vector of shape (m,).

• rank (int, optional) – Truncation rank for regularized TLS. If None, uses full rank. Useful when A is rank-deficient.

Returns:

result – Named tuple containing: - x: Estimated parameters of shape (n,) - residuals_A: Corrections to A of shape (m, n) - residuals_b: Corrections to b of shape (m,) - rank: Effective rank used

Return type:

TLSResult

Examples

>>> A = np.array([[1, 1], [2, 1], [3, 1]])
>>> b = np.array([2.1, 2.9, 4.1])  # Noisy measurements
>>> result = total_least_squares(A, b)

Notes

TLS is appropriate when both the independent and dependent variables are subject to measurement error. It minimizes the orthogonal distance from the data points to the fitted model.

The solution is computed using the SVD of the augmented matrix [A | b].

References

[1]

S. Van Huffel and J. Vandewalle, “The Total Least Squares Problem: Computational Aspects and Analysis,” SIAM, 1991.

pytcl.static_estimation.generalized_least_squares(A, b, Sigma)

Generalized least squares (GLS) estimation.

Solves the linear model with correlated errors:

b = A @ x + e, where E[e e^T] = Sigma

This is equivalent to weighted least squares with W = Sigma^{-1}.

Parameters:

• A (array_like) – Design matrix of shape (m, n).

• b (array_like) – Observation vector of shape (m,).

• Sigma (array_like) – Error covariance matrix of shape (m, m).

Returns:

result – Same as weighted_least_squares result.

Return type:

WLSResult

Examples

>>> A = np.array([[1, 1], [1, 2], [1, 3]])
>>> b = np.array([1, 2, 2])
>>> Sigma = np.array([[1, 0.5, 0], [0.5, 1, 0.5], [0, 0.5, 1]])
>>> result = generalized_least_squares(A, b, Sigma)

Notes

GLS is the BLUE (Best Linear Unbiased Estimator) when the error covariance structure is known.

pytcl.static_estimation.recursive_least_squares(x_prev, P_prev, a, y, forgetting_factor=1.0)

Recursive least squares (RLS) update.

Updates parameter estimates when a new observation arrives, without reprocessing all previous data.

Parameters:

• x_prev (array_like) – Previous parameter estimate of shape (n,).

• P_prev (array_like) – Previous covariance matrix of shape (n, n).

• a (array_like) – New regressor vector of shape (n,).

• y (float) – New observation.

• forgetting_factor (float, optional) – Forgetting factor in (0, 1]. Values < 1 give more weight to recent observations. Default is 1.0 (no forgetting).

Returns:

• x (ndarray) – Updated parameter estimate.

• P (ndarray) – Updated covariance matrix.

Return type:

tuple[ndarray[tuple[Any, …], dtype[floating]], ndarray[tuple[Any, …], dtype[floating]]]

Examples

>>> x = np.zeros(2)  # Initial estimate
>>> P = np.eye(2) * 100  # High initial uncertainty
>>> # Process observations one at a time
>>> x, P = recursive_least_squares(x, P, np.array([1, 1]), 2.0)
>>> x, P = recursive_least_squares(x, P, np.array([1, 2]), 3.0)

Notes

RLS is equivalent to the Kalman filter for the static parameter estimation problem. The forgetting factor introduces exponential weighting of past data.

pytcl.static_estimation.ridge_regression(A, b, alpha=1.0)

Ridge regression (L2-regularized least squares).

Solves the regularized problem:

min_x ||A @ x - b||_2^2 + alpha * ||x||_2^2

Parameters:

• A (array_like) – Design matrix of shape (m, n).

• b (array_like) – Observation vector of shape (m,).

• alpha (float, optional) – Regularization parameter. Larger values give more regularization. Default is 1.0.

Returns:

x – Regularized parameter estimate.

Return type:

ndarray

Examples

>>> A = np.array([[1, 1], [1, 2], [1, 3]])
>>> b = np.array([1, 2, 2])
>>> x = ridge_regression(A, b, alpha=0.1)

Notes

Ridge regression shrinks parameters toward zero and is useful when: - The design matrix is ill-conditioned - There are more parameters than observations (n > m) - Regularization is desired to prevent overfitting

class pytcl.static_estimation.RobustResult(x, residuals, weights, scale, n_iter, converged)

Bases: NamedTuple

Result of robust estimation.

x

Estimated parameters.

Type:

ndarray

residuals

Residual vector.

Type:

ndarray

weights

Final weights for each observation.

Type:

ndarray

scale

Estimated scale of residuals.

Type:

float

n_iter

Number of iterations performed.

Type:

int

converged

Whether the algorithm converged.

Type:

bool

x: ndarray[tuple[Any, ...], dtype[floating]]

Alias for field number 0

residuals: ndarray[tuple[Any, ...], dtype[floating]]

Alias for field number 1

weights: ndarray[tuple[Any, ...], dtype[floating]]

Alias for field number 2

scale: float

Alias for field number 3

n_iter: int

Alias for field number 4

converged: bool

Alias for field number 5

class pytcl.static_estimation.RANSACResult(x, inliers, n_inliers, residuals, n_iter, best_score)

Bases: NamedTuple

Result of RANSAC estimation.

x

Estimated parameters from best model.

Type:

ndarray

inliers

Boolean mask of inlier points.

Type:

ndarray

n_inliers

Number of inliers.

Type:

int

residuals

Residuals for all points.

Type:

ndarray

n_iter

Number of iterations performed.

Type:

int

best_score

Score of the best model found.

Type:

float

x: ndarray[tuple[Any, ...], dtype[floating]]

Alias for field number 0

inliers: ndarray[tuple[Any, ...], dtype[bool]]

Alias for field number 1

n_inliers: int

Alias for field number 2

residuals: ndarray[tuple[Any, ...], dtype[floating]]

Alias for field number 3

n_iter: int

Alias for field number 4

best_score: float

Alias for field number 5

pytcl.static_estimation.huber_weight(r, c=1.345)

Huber weight function.

Parameters:

• r (array_like) – Standardized residuals.

• c (float, optional) – Tuning constant. Default 1.345 gives 95% efficiency for normal distribution.

Returns:

weights – Weights for each residual.

Return type:

ndarray

Notes

The Huber weight function is:

w(r) = 1 if |r| <= c w(r) = c / |r| if |r| > c

Examples

>>> r = np.array([0.5, 1.0, 2.0, 5.0])  # Standardized residuals
>>> w = huber_weight(r, c=1.345)
>>> w[0]  # Small residual gets weight 1
1.0
>>> w[3] < 0.5  # Large residual gets reduced weight
True

pytcl.static_estimation.huber_rho(r, c=1.345)

Huber rho (loss) function.

Parameters:

• r (array_like) – Standardized residuals.

• c (float, optional) – Tuning constant.

Returns:

rho – Loss values.

Return type:

ndarray

Notes

The Huber rho function is:

rho(r) = r^2 / 2 if |r| <= c rho(r) = c * |r| - c^2/2 if |r| > c

Examples

>>> r = np.array([0.5, 1.0, 2.0])
>>> rho = huber_rho(r, c=1.345)
>>> rho[0]  # Small residual: r^2/2
0.125

pytcl.static_estimation.tukey_weight(r, c=4.685)

Tukey bisquare weight function.

Parameters:

• r (array_like) – Standardized residuals.

• c (float, optional) – Tuning constant. Default 4.685 gives 95% efficiency for normal distribution.

Returns:

weights – Weights for each residual.

Return type:

ndarray

Notes

The Tukey bisquare weight function is:

w(r) = (1 - (r/c)^2)^2 if |r| <= c w(r) = 0 if |r| > c

This provides complete rejection of large outliers.

Examples

>>> r = np.array([0.5, 2.0, 5.0, 10.0])
>>> w = tukey_weight(r, c=4.685)
>>> w[0] > 0.9  # Small residual gets high weight
True
>>> w[3]  # Large residual completely rejected
0.0

pytcl.static_estimation.tukey_rho(r, c=4.685)

Tukey bisquare rho (loss) function.

Parameters:

• r (array_like) – Standardized residuals.

• c (float, optional) – Tuning constant.

Returns:

rho – Loss values.

Return type:

ndarray

Notes

The Tukey rho function is:

rho(r) = c^2/6 * (1 - (1 - (r/c)^2)^3) if |r| <= c rho(r) = c^2/6 if |r| > c

Examples

>>> r = np.array([0.0, 2.0, 10.0])
>>> rho = tukey_rho(r, c=4.685)
>>> rho[0]  # Zero residual
0.0
>>> rho[2] == rho[2]  # Large residuals saturate at c^2/6
True

pytcl.static_estimation.cauchy_weight(r, c=2.385)

Cauchy weight function.

Parameters:

• r (array_like) – Standardized residuals.

• c (float, optional) – Tuning constant.

Returns:

weights – Weights for each residual.

Return type:

ndarray

Notes

The Cauchy weight function is:

w(r) = 1 / (1 + (r/c)^2)

Examples

>>> r = np.array([0.0, 1.0, 5.0])
>>> w = cauchy_weight(r, c=2.385)
>>> w[0]  # Zero residual gets weight 1
1.0
>>> 0 < w[2] < 1  # Large residuals get reduced weight (but never zero)
True

pytcl.static_estimation.mad(residuals, c=1.4826)

Median Absolute Deviation (MAD) scale estimator.

Parameters:

• residuals (array_like) – Residual vector.

• c (float, optional) – Consistency constant. Default 1.4826 makes MAD consistent for normal distribution.

Returns:

scale – Estimated scale.

Return type:

float

Notes

MAD = c * median(|r - median(r)|)

This is a robust scale estimator with 50% breakdown point.

Examples

>>> residuals = np.array([1.0, 1.1, 0.9, 1.0, 100.0])  # One outlier
>>> scale = mad(residuals)
>>> scale < 1.0  # Robust to the outlier
True

pytcl.static_estimation.tau_scale(residuals, c1=4.5, c2=3.0)

Tau scale estimator.

Parameters:

• residuals (array_like) – Residual vector.

• c1 (float) – Tuning constants.

• c2 (float) – Tuning constants.

Returns:

scale – Estimated scale.

Return type:

float

Notes

Tau scale combines high breakdown point with efficiency.

Examples

>>> residuals = np.array([1.0, 1.1, 0.9, 1.0, 1.2, 100.0])  # One outlier
>>> scale = tau_scale(residuals)
>>> scale < 10.0  # Robust to the outlier
True

pytcl.static_estimation.irls(A, b, weight_func=<function huber_weight>, scale_func=<function mad>, max_iter=50, tol=1e-06)

Iteratively Reweighted Least Squares (IRLS).

General M-estimator using IRLS algorithm.

Parameters:

• A (array_like) – Design matrix of shape (m, n).

• b (array_like) – Observation vector of shape (m,).

• weight_func (callable, optional) – Weight function w(r) for standardized residuals. Default is Huber weight.

• scale_func (callable, optional) – Scale estimation function. Default is MAD.

• max_iter (int, optional) – Maximum number of iterations. Default 50.

• tol (float, optional) – Convergence tolerance for parameter change. Default 1e-6.

Returns:

result – Robust estimation result.

Return type:

RobustResult

Examples

>>> A = np.array([[1, 1], [1, 2], [1, 3], [1, 4], [1, 10]])
>>> b = np.array([2, 3, 4, 5, 100])  # Last point is outlier
>>> result = irls(A, b)
>>> result.x  # Should be close to [1, 1] ignoring outlier

Notes

IRLS iteratively solves weighted least squares problems, updating weights based on the current residuals.

pytcl.static_estimation.huber_regression(A, b, c=1.345, max_iter=50, tol=1e-06)

Huber robust regression.

M-estimation using Huber’s weight function.

Parameters:

• A (array_like) – Design matrix of shape (m, n).

• b (array_like) – Observation vector of shape (m,).

• c (float, optional) – Tuning constant. Default 1.345 gives 95% efficiency.

• max_iter (int, optional) – Maximum iterations. Default 50.

• tol (float, optional) – Convergence tolerance. Default 1e-6.

Returns:

result – Robust estimation result.

Return type:

RobustResult

Examples

>>> A = np.array([[1, 1], [1, 2], [1, 3]])
>>> b = np.array([2.1, 2.9, 100])  # Outlier in last observation
>>> result = huber_regression(A, b)

Notes

Huber regression provides a balance between efficiency for Gaussian errors and resistance to outliers.

pytcl.static_estimation.tukey_regression(A, b, c=4.685, max_iter=50, tol=1e-06)

Tukey bisquare robust regression.

M-estimation using Tukey’s bisquare weight function.

Parameters:

• A (array_like) – Design matrix of shape (m, n).

• b (array_like) – Observation vector of shape (m,).

• c (float, optional) – Tuning constant. Default 4.685 gives 95% efficiency.

• max_iter (int, optional) – Maximum iterations. Default 50.

• tol (float, optional) – Convergence tolerance. Default 1e-6.

Returns:

result – Robust estimation result.

Return type:

RobustResult

Examples

>>> A = np.array([[1, 1], [1, 2], [1, 3]])
>>> b = np.array([2.1, 2.9, 100])  # Outlier in last observation
>>> result = tukey_regression(A, b)

Notes

Tukey bisquare provides complete rejection of large outliers (zero weight for residuals > c), making it more robust than Huber for gross outliers.

pytcl.static_estimation.ransac(A, b, min_samples=None, residual_threshold=None, max_trials=100, stop_n_inliers=None, stop_score=None, random_state=None)

RANdom SAmple Consensus (RANSAC) regression.

Robust regression by iteratively selecting random subsets, fitting models, and identifying inliers.

Parameters:

• A (array_like) – Design matrix of shape (m, n).

• b (array_like) – Observation vector of shape (m,).

• min_samples (int, optional) – Minimum number of samples for fitting. Default is n (number of features).

• residual_threshold (float, optional) – Maximum residual for a point to be considered inlier. Default is MAD of initial OLS residuals.

• max_trials (int, optional) – Maximum number of random trials. Default 100.

• stop_n_inliers (int, optional) – Stop early if this many inliers found. Default None.

• stop_score (float, optional) – Stop early if score exceeds this. Default None.

• random_state (int, optional) – Random seed for reproducibility.

Returns:

result – RANSAC result with best model and inlier mask.

Return type:

RANSACResult

Examples

>>> rng = np.random.default_rng(42)
>>> A = np.column_stack([np.ones(100), np.arange(100)])
>>> b = 2 + 3 * np.arange(100) + rng.normal(0, 1, 100)
>>> # Add outliers
>>> b[90:] = 1000
>>> result = ransac(A, b, random_state=42)
>>> result.x  # Should be close to [2, 3]

Notes

RANSAC is particularly effective when the fraction of outliers is large (>25%), where M-estimators may struggle.

The algorithm: 1. Randomly select min_samples points 2. Fit model to selected points 3. Count inliers (points with residual < threshold) 4. If best model found, save it 5. Repeat for max_trials 6. Refit model using all inliers of best model

References

[1]

M. A. Fischler and R. C. Bolles, “Random Sample Consensus,” Communications of the ACM, 1981.

pytcl.static_estimation.ransac_n_trials(n_samples, n_outliers, min_samples, probability=0.99)

Compute number of RANSAC trials needed.

Parameters:

• n_samples (int) – Total number of samples.

• n_outliers (int) – Expected number of outliers.

• min_samples (int) – Number of samples per trial.

• probability (float, optional) – Desired probability of success. Default 0.99.

Returns:

n_trials – Number of trials needed.

Return type:

int

Examples

>>> ransac_n_trials(100, 30, 2)  # 30% outliers, 2 samples per trial
10

Notes

Formula: k = log(1 - p) / log(1 - (1 - e)^n) where: - p = probability of success - e = outlier ratio - n = min_samples

class pytcl.static_estimation.MLResult(theta, log_likelihood, fisher_info, covariance, n_iter, converged)

Bases: NamedTuple

Result of maximum likelihood estimation.

theta

Estimated parameters.

Type:

ndarray

log_likelihood

Log-likelihood at the estimate.

Type:

float

fisher_info

Fisher information matrix at the estimate.

Type:

ndarray

covariance

Estimated covariance (inverse Fisher information).

Type:

ndarray

n_iter

Number of iterations (for iterative methods).

Type:

int

converged

Whether the optimization converged.

Type:

bool

theta: ndarray[tuple[Any, ...], dtype[floating]]

Alias for field number 0

log_likelihood: float

Alias for field number 1

fisher_info: ndarray[tuple[Any, ...], dtype[floating]]

Alias for field number 2

covariance: ndarray[tuple[Any, ...], dtype[floating]]

Alias for field number 3

n_iter: int

Alias for field number 4

converged: bool

Alias for field number 5

class pytcl.static_estimation.CRBResult(crb_matrix, variances, std_bounds, fisher_info, is_efficient)

Bases: NamedTuple

Result of Cramer-Rao bound computation.

crb_matrix

Cramer-Rao bound matrix (inverse Fisher information).

Type:

ndarray

variances

Diagonal elements (variance bounds for each parameter).

Type:

ndarray

std_bounds

Square root of variances (standard deviation bounds).

Type:

ndarray

fisher_info

Fisher information matrix used.

Type:

ndarray

is_efficient

Boolean mask indicating which estimators achieve the bound.

Type:

ndarray or None

crb_matrix: ndarray[tuple[Any, ...], dtype[floating]]

Alias for field number 0

variances: ndarray[tuple[Any, ...], dtype[floating]]

Alias for field number 1

std_bounds: ndarray[tuple[Any, ...], dtype[floating]]

Alias for field number 2

fisher_info: ndarray[tuple[Any, ...], dtype[floating]]

Alias for field number 3

is_efficient: ndarray[tuple[Any, ...], dtype[bool]] | None

Alias for field number 4

pytcl.static_estimation.fisher_information_numerical(log_likelihood, theta, h=1e-05)

Compute Fisher information matrix numerically.

Uses the negative expected Hessian of the log-likelihood.

Parameters:

• log_likelihood (callable) – Log-likelihood function L(theta).

• theta (array_like) – Parameter vector at which to evaluate.

• h (float, optional) – Step size for numerical differentiation. Default 1e-5.

Returns:

fisher – Fisher information matrix of shape (n_params, n_params).

Return type:

ndarray

Examples

>>> def log_lik(theta):
...     return -0.5 * np.sum((data - theta[0])**2 / theta[1])
>>> theta = np.array([0.0, 1.0])
>>> F = fisher_information_numerical(log_lik, theta)

Notes

For a single observation, the Fisher information is:

I(theta) = -E[d^2 log p(x|theta) / d theta^2]

This function approximates the Hessian using central differences.

pytcl.static_estimation.fisher_information_gaussian(jacobian, noise_cov)

Fisher information for Gaussian measurement model.

For the model y = h(theta) + noise, where noise ~ N(0, R), the Fisher information is J^T R^{-1} J.

Parameters:

• jacobian (array_like) – Jacobian matrix dh/dtheta of shape (m, n).

• noise_cov (array_like) – Measurement noise covariance R of shape (m, m).

Returns:

fisher – Fisher information matrix of shape (n, n).

Return type:

ndarray

Examples

>>> H = np.array([[1, 0], [0, 1], [1, 1]])  # 3 measurements, 2 params
>>> R = np.eye(3) * 0.1
>>> F = fisher_information_gaussian(H, R)

Notes

This is the standard result for linear Gaussian models and provides the information content of measurements about parameters.

pytcl.static_estimation.fisher_information_exponential_family(sufficient_stats, theta, data, h=1e-05)

Fisher information for exponential family distributions.

For exponential family: p(x|theta) = h(x) exp(eta(theta)^T T(x) - A(theta)) The Fisher information equals the covariance of sufficient statistics.

Parameters:

• sufficient_stats (callable) – Function T(x, theta) returning sufficient statistics.

• theta (array_like) – Natural parameters.

• data (array_like) – Observed data of shape (n_samples, …).

• h (float, optional) – Step size for numerical differentiation.

Returns:

fisher – Fisher information matrix.

Return type:

ndarray

Notes

For exponential families, I(theta) = Var[T(X)] = d^2 A(theta) / d theta^2, where A(theta) is the log-partition function.

Examples

>>> def suff_stats(x, theta):
...     return np.array([x, x**2])  # Mean and second moment
>>> data = np.random.normal(0, 1, 100)
>>> theta = np.array([0.0, 1.0])
>>> F = fisher_information_exponential_family(suff_stats, theta, data)

pytcl.static_estimation.observed_fisher_information(log_likelihood, theta, h=1e-05)

Compute observed Fisher information (negative Hessian).

Unlike the expected Fisher information, this uses the actual observed data rather than the expectation.

Parameters:

• log_likelihood (callable) – Log-likelihood function.

• theta (array_like) – Parameter estimate.

• h (float, optional) – Step size for numerical differentiation.

Returns:

observed_fisher – Observed Fisher information matrix.

Return type:

ndarray

Notes

The observed Fisher information is often more accurate for finite samples and is asymptotically equivalent to the expected Fisher information.

Examples

>>> data = np.array([1.2, 0.8, 1.1, 0.9, 1.0])
>>> def log_lik(theta):
...     return -0.5 * np.sum((data - theta[0])**2 / theta[1]) - 2.5 * np.log(theta[1])
>>> theta = np.array([1.0, 0.1])
>>> F_obs = observed_fisher_information(log_lik, theta)

pytcl.static_estimation.cramer_rao_bound(fisher_info, estimator_variance=None)

Compute Cramer-Rao lower bound on estimator variance.

Parameters:

• fisher_info (array_like) – Fisher information matrix of shape (n, n).

• estimator_variance (array_like, optional) – Actual estimator variance for efficiency check.

Returns:

result – Named tuple with CRB matrix, variances, and efficiency info.

Return type:

CRBResult

Examples

>>> F = np.array([[10, 0], [0, 5]])  # Fisher info
>>> result = cramer_rao_bound(F)
>>> result.variances  # Minimum achievable variances
array([0.1, 0.2])

Notes

The Cramer-Rao bound states that for any unbiased estimator:

Var(theta_hat) >= I(theta)^{-1}

The bound is achieved by efficient estimators, such as the MLE for exponential family distributions.

pytcl.static_estimation.cramer_rao_bound_biased(fisher_info, bias_gradient)

Cramer-Rao bound for biased estimators.

Parameters:

• fisher_info (array_like) – Fisher information matrix.

• bias_gradient (array_like) – Gradient of the bias with respect to theta.

Returns:

crb – Cramer-Rao bound matrix for the biased estimator.

Return type:

ndarray

Notes

For a biased estimator with bias b(theta), the CRB becomes:

Var(theta_hat) >= (I + db/dtheta) I^{-1} (I + db/dtheta)^T

Examples

>>> F = np.array([[10.0, 0], [0, 5.0]])  # Fisher info
>>> db = np.array([[0.1, 0], [0, 0.2]])  # Bias gradient
>>> crb_biased = cramer_rao_bound_biased(F, db)
>>> crb_biased.shape
(2, 2)

pytcl.static_estimation.efficiency(estimator_variance, crb)

Compute estimator efficiency relative to CRB.

Parameters:

• estimator_variance (array_like) – Variance of the estimator (scalar or diagonal).

• crb (array_like) – Cramer-Rao bound (scalar or diagonal).

Returns:

eff – Efficiency values in [0, 1]. Value of 1 means efficient.

Return type:

ndarray

Examples

>>> var_est = np.array([0.12, 0.25])
>>> crb = np.array([0.1, 0.2])
>>> efficiency(var_est, crb)
array([0.833, 0.8])

pytcl.static_estimation.mle_newton_raphson(log_likelihood, score, theta_init, hessian=None, max_iter=100, tol=1e-08, h=1e-05)

Maximum likelihood estimation using Newton-Raphson.

Parameters:

• log_likelihood (callable) – Log-likelihood function L(theta).

• score (callable) – Score function (gradient of log-likelihood).

• theta_init (array_like) – Initial parameter guess.

• hessian (callable, optional) – Hessian function. If None, computed numerically.

• max_iter (int, optional) – Maximum iterations. Default 100.

• tol (float, optional) – Convergence tolerance. Default 1e-8.

• h (float, optional) – Step size for numerical Hessian.

Returns:

result – MLE result with estimate, Fisher info, and covariance.

Return type:

MLResult

Examples

>>> def log_lik(theta):
...     return -0.5 * np.sum((data - theta[0])**2)
>>> def score(theta):
...     return np.array([np.sum(data - theta[0])])
>>> result = mle_newton_raphson(log_lik, score, np.array([0.0]))

Notes

Newton-Raphson update: theta_{n+1} = theta_n - H^{-1} @ score where H is the Hessian of the log-likelihood.

pytcl.static_estimation.mle_scoring(log_likelihood, score, fisher_info, theta_init, max_iter=100, tol=1e-08)

Maximum likelihood estimation using Fisher scoring.

Uses expected Fisher information instead of observed Hessian, which can be more stable.

Parameters:

• log_likelihood (callable) – Log-likelihood function.

• score (callable) – Score function (gradient).

• fisher_info (callable) – Fisher information function.

• theta_init (array_like) – Initial parameter guess.

• max_iter (int, optional) – Maximum iterations.

• tol (float, optional) – Convergence tolerance.

Returns:

result – MLE result.

Return type:

MLResult

Notes

Fisher scoring update: theta_{n+1} = theta_n + I(theta_n)^{-1} @ score This is equivalent to Newton-Raphson when I(theta) = -E[H].

Examples

>>> data = np.array([1.0, 1.1, 0.9, 1.2, 0.8])
>>> def log_lik(theta):
...     return -0.5 * len(data) * np.log(2*np.pi) - np.sum((data - theta[0])**2) / 2
>>> def score(theta):
...     return np.array([np.sum(data - theta[0])])
>>> def fisher(theta):
...     return np.array([[len(data)]])
>>> result = mle_scoring(log_lik, score, fisher, np.array([0.0]))

pytcl.static_estimation.mle_gaussian(data, estimate_mean=True, estimate_variance=True)

Closed-form MLE for Gaussian distribution.

Parameters:

• data (array_like) – Observed data of shape (n_samples,) or (n_samples, n_features).

• estimate_mean (bool, optional) – Whether to estimate mean. Default True.

• estimate_variance (bool, optional) – Whether to estimate variance. Default True.

Returns:

result – MLE result with mean and/or variance estimates.

Return type:

MLResult

Examples

>>> data = np.random.normal(5, 2, 1000)
>>> result = mle_gaussian(data)
>>> result.theta  # [mean, variance]

pytcl.static_estimation.aic(log_likelihood, n_params)

Akaike Information Criterion.

Parameters:

• log_likelihood (float) – Log-likelihood at the MLE.

• n_params (int) – Number of parameters.

Returns:

aic – AIC value (lower is better).

Return type:

float

Notes

AIC = -2 * log_likelihood + 2 * n_params

Examples

>>> log_lik = -100.0
>>> n_params = 3
>>> aic(log_lik, n_params)
206.0

pytcl.static_estimation.bic(log_likelihood, n_params, n_samples)

Bayesian Information Criterion.

Parameters:

• log_likelihood (float) – Log-likelihood at the MLE.

• n_params (int) – Number of parameters.

• n_samples (int) – Number of samples.

Returns:

bic – BIC value (lower is better).

Return type:

float

Notes

BIC = -2 * log_likelihood + n_params * log(n_samples)

Examples

>>> log_lik = -100.0
>>> bic(log_lik, n_params=3, n_samples=100)
213.81551055796427

pytcl.static_estimation.aicc(log_likelihood, n_params, n_samples)

Corrected Akaike Information Criterion.

Parameters:

• log_likelihood (float) – Log-likelihood at the MLE.

• n_params (int) – Number of parameters.

• n_samples (int) – Number of samples.

Returns:

aicc – AICc value (lower is better).

Return type:

float

Notes

AICc adds a correction for small sample sizes: AICc = AIC + 2*k*(k+1)/(n-k-1)

Examples

>>> log_lik = -50.0
>>> aicc(log_lik, n_params=3, n_samples=20)
109.5

Least Squares

Ordinary, weighted, total, generalized, and recursive least squares.

Least squares estimation methods.

This module provides ordinary, weighted, and total least squares estimators commonly used in tracking for state estimation and model fitting.

References

[1]

S. Van Huffel and J. Vandewalle, “The Total Least Squares Problem: Computational Aspects and Analysis,” SIAM, 1991.

[2]

G. H. Golub and C. F. Van Loan, “Matrix Computations,” 4th ed., Johns Hopkins University Press, 2013.

class pytcl.static_estimation.least_squares.LSResult(x, residuals, rank, singular_values)

Bases: NamedTuple

Result of least squares estimation.

x

Estimated parameters.

Type:

ndarray

residuals

Residual vector (y - A @ x).

Type:

ndarray

rank

Effective rank of the design matrix.

Type:

int

singular_values

Singular values of the design matrix.

Type:

ndarray

x: ndarray[tuple[Any, ...], dtype[floating]]

Alias for field number 0

residuals: ndarray[tuple[Any, ...], dtype[floating]]

Alias for field number 1

rank: int

Alias for field number 2

singular_values: ndarray[tuple[Any, ...], dtype[floating]]

Alias for field number 3

class pytcl.static_estimation.least_squares.WLSResult(x, residuals, covariance, weighted_residual_sum)

Bases: NamedTuple

Result of weighted least squares estimation.

x

Estimated parameters.

Type:

ndarray

residuals

Residual vector.

Type:

ndarray

covariance

Estimated covariance of parameters.

Type:

ndarray

weighted_residual_sum

Sum of weighted squared residuals.

Type:

float

x: ndarray[tuple[Any, ...], dtype[floating]]

Alias for field number 0

residuals: ndarray[tuple[Any, ...], dtype[floating]]

Alias for field number 1

covariance: ndarray[tuple[Any, ...], dtype[floating]]

Alias for field number 2

weighted_residual_sum: float

Alias for field number 3

class pytcl.static_estimation.least_squares.TLSResult(x, residuals_A, residuals_b, rank)

Bases: NamedTuple

Result of total least squares estimation.

x

Estimated parameters.

Type:

ndarray

residuals_A

Corrections to the design matrix A.

Type:

ndarray

residuals_b

Corrections to the observation vector b.

Type:

ndarray

rank

Effective rank used in the solution.

Type:

int

x: ndarray[tuple[Any, ...], dtype[floating]]

Alias for field number 0

residuals_A: ndarray[tuple[Any, ...], dtype[floating]]

Alias for field number 1

residuals_b: ndarray[tuple[Any, ...], dtype[floating]]

Alias for field number 2

rank: int

Alias for field number 3

pytcl.static_estimation.least_squares.ordinary_least_squares(A, b, rcond=None)

Ordinary least squares estimation.

Solves the linear least squares problem:

min_x ||A @ x - b||_2

Parameters:

• A (array_like) – Design matrix of shape (m, n) where m >= n.

• b (array_like) – Observation vector of shape (m,) or (m, k).

• rcond (float, optional) – Cutoff for small singular values. Values smaller than rcond * largest_singular_value are treated as zero. Default is machine precision * max(m, n).

Returns:

result – Named tuple containing: - x: Estimated parameters of shape (n,) or (n, k) - residuals: Residual vector - rank: Effective rank of A - singular_values: Singular values of A

Return type:

LSResult

Examples

>>> A = np.array([[1, 1], [1, 2], [1, 3]])
>>> b = np.array([1, 2, 2])
>>> result = ordinary_least_squares(A, b)
>>> result.x  # Fitted line parameters
array([0.66666667, 0.5       ])

Notes

Uses SVD-based solution for numerical stability.

pytcl.static_estimation.least_squares.weighted_least_squares(A, b, W=None, weights=None)

Weighted least squares estimation.

Solves the weighted linear least squares problem:

min_x (b - A @ x)^T W (b - A @ x)

Parameters:

• A (array_like) – Design matrix of shape (m, n).

• b (array_like) – Observation vector of shape (m,).

• W (array_like, optional) – Weight matrix of shape (m, m). If provided, weights is ignored. Should be positive definite.

• weights (array_like, optional) – Diagonal weights of shape (m,). Used only if W is None. Equivalent to W = diag(weights).

Returns:

result – Named tuple containing: - x: Estimated parameters of shape (n,) - residuals: Residual vector - covariance: Estimated parameter covariance - weighted_residual_sum: Weighted sum of squared residuals

Return type:

WLSResult

Examples

>>> A = np.array([[1, 1], [1, 2], [1, 3]])
>>> b = np.array([1, 2, 2])
>>> weights = np.array([1, 2, 1])  # Higher weight on middle observation
>>> result = weighted_least_squares(A, b, weights=weights)

Notes

The covariance of the estimated parameters is (A^T W A)^{-1}.

For measurement noise with covariance R, use W = R^{-1} to get the minimum variance unbiased estimator (MVUE).

pytcl.static_estimation.least_squares.total_least_squares(A, b, rank=None)

Total least squares (TLS) estimation.

Solves the errors-in-variables problem where both the design matrix and observations may have errors:

min_{E, r} ||[E | r]||_F subject to (A + E) @ x = b + r

Parameters:

• A (array_like) – Design matrix of shape (m, n) where m >= n.

• b (array_like) – Observation vector of shape (m,).

• rank (int, optional) – Truncation rank for regularized TLS. If None, uses full rank. Useful when A is rank-deficient.

Returns:

result – Named tuple containing: - x: Estimated parameters of shape (n,) - residuals_A: Corrections to A of shape (m, n) - residuals_b: Corrections to b of shape (m,) - rank: Effective rank used

Return type:

TLSResult

Examples

>>> A = np.array([[1, 1], [2, 1], [3, 1]])
>>> b = np.array([2.1, 2.9, 4.1])  # Noisy measurements
>>> result = total_least_squares(A, b)

Notes

TLS is appropriate when both the independent and dependent variables are subject to measurement error. It minimizes the orthogonal distance from the data points to the fitted model.

The solution is computed using the SVD of the augmented matrix [A | b].

References

[1]

S. Van Huffel and J. Vandewalle, “The Total Least Squares Problem: Computational Aspects and Analysis,” SIAM, 1991.

pytcl.static_estimation.least_squares.generalized_least_squares(A, b, Sigma)

Generalized least squares (GLS) estimation.

Solves the linear model with correlated errors:

b = A @ x + e, where E[e e^T] = Sigma

This is equivalent to weighted least squares with W = Sigma^{-1}.

Parameters:

• A (array_like) – Design matrix of shape (m, n).

• b (array_like) – Observation vector of shape (m,).

• Sigma (array_like) – Error covariance matrix of shape (m, m).

Returns:

result – Same as weighted_least_squares result.

Return type:

WLSResult

Examples

>>> A = np.array([[1, 1], [1, 2], [1, 3]])
>>> b = np.array([1, 2, 2])
>>> Sigma = np.array([[1, 0.5, 0], [0.5, 1, 0.5], [0, 0.5, 1]])
>>> result = generalized_least_squares(A, b, Sigma)

Notes

GLS is the BLUE (Best Linear Unbiased Estimator) when the error covariance structure is known.

pytcl.static_estimation.least_squares.recursive_least_squares(x_prev, P_prev, a, y, forgetting_factor=1.0)

Recursive least squares (RLS) update.

Updates parameter estimates when a new observation arrives, without reprocessing all previous data.

Parameters:

• x_prev (array_like) – Previous parameter estimate of shape (n,).

• P_prev (array_like) – Previous covariance matrix of shape (n, n).

• a (array_like) – New regressor vector of shape (n,).

• y (float) – New observation.

• forgetting_factor (float, optional) – Forgetting factor in (0, 1]. Values < 1 give more weight to recent observations. Default is 1.0 (no forgetting).

Returns:

• x (ndarray) – Updated parameter estimate.

• P (ndarray) – Updated covariance matrix.

Return type:

tuple[ndarray[tuple[Any, …], dtype[floating]], ndarray[tuple[Any, …], dtype[floating]]]

Examples

>>> x = np.zeros(2)  # Initial estimate
>>> P = np.eye(2) * 100  # High initial uncertainty
>>> # Process observations one at a time
>>> x, P = recursive_least_squares(x, P, np.array([1, 1]), 2.0)
>>> x, P = recursive_least_squares(x, P, np.array([1, 2]), 3.0)

Notes

RLS is equivalent to the Kalman filter for the static parameter estimation problem. The forgetting factor introduces exponential weighting of past data.

pytcl.static_estimation.least_squares.ridge_regression(A, b, alpha=1.0)

Ridge regression (L2-regularized least squares).

Solves the regularized problem:

min_x ||A @ x - b||_2^2 + alpha * ||x||_2^2

Parameters:

• A (array_like) – Design matrix of shape (m, n).

• b (array_like) – Observation vector of shape (m,).

• alpha (float, optional) – Regularization parameter. Larger values give more regularization. Default is 1.0.

Returns:

x – Regularized parameter estimate.

Return type:

ndarray

Examples

>>> A = np.array([[1, 1], [1, 2], [1, 3]])
>>> b = np.array([1, 2, 2])
>>> x = ridge_regression(A, b, alpha=0.1)

Notes

Ridge regression shrinks parameters toward zero and is useful when: - The design matrix is ill-conditioned - There are more parameters than observations (n > m) - Regularization is desired to prevent overfitting

Maximum Likelihood

Maximum likelihood estimation and Fisher information.

Maximum Likelihood Estimation and Information Theory.

This module provides tools for maximum likelihood estimation, Fisher information computation, and Cramer-Rao bound analysis.

References

[1]

S. M. Kay, “Fundamentals of Statistical Signal Processing: Estimation Theory,” Prentice Hall, 1993.

[2]

H. L. Van Trees, “Detection, Estimation, and Modulation Theory,” Wiley, 2001.

class pytcl.static_estimation.maximum_likelihood.MLResult(theta, log_likelihood, fisher_info, covariance, n_iter, converged)

Bases: NamedTuple

Result of maximum likelihood estimation.

theta

Estimated parameters.

Type:

ndarray

log_likelihood

Log-likelihood at the estimate.

Type:

float

fisher_info

Fisher information matrix at the estimate.

Type:

ndarray

covariance

Estimated covariance (inverse Fisher information).

Type:

ndarray

n_iter

Number of iterations (for iterative methods).

Type:

int

converged

Whether the optimization converged.

Type:

bool

theta: ndarray[tuple[Any, ...], dtype[floating]]

Alias for field number 0

log_likelihood: float

Alias for field number 1

fisher_info: ndarray[tuple[Any, ...], dtype[floating]]

Alias for field number 2

covariance: ndarray[tuple[Any, ...], dtype[floating]]

Alias for field number 3

n_iter: int

Alias for field number 4

converged: bool

Alias for field number 5

class pytcl.static_estimation.maximum_likelihood.CRBResult(crb_matrix, variances, std_bounds, fisher_info, is_efficient)

Bases: NamedTuple

Result of Cramer-Rao bound computation.

crb_matrix

Cramer-Rao bound matrix (inverse Fisher information).

Type:

ndarray

variances

Diagonal elements (variance bounds for each parameter).

Type:

ndarray

std_bounds

Square root of variances (standard deviation bounds).

Type:

ndarray

fisher_info

Fisher information matrix used.

Type:

ndarray

is_efficient

Boolean mask indicating which estimators achieve the bound.

Type:

ndarray or None

crb_matrix: ndarray[tuple[Any, ...], dtype[floating]]

Alias for field number 0

variances: ndarray[tuple[Any, ...], dtype[floating]]

Alias for field number 1

std_bounds: ndarray[tuple[Any, ...], dtype[floating]]

Alias for field number 2

fisher_info: ndarray[tuple[Any, ...], dtype[floating]]

Alias for field number 3

is_efficient: ndarray[tuple[Any, ...], dtype[bool]] | None

Alias for field number 4

pytcl.static_estimation.maximum_likelihood.fisher_information_numerical(log_likelihood, theta, h=1e-05)

Compute Fisher information matrix numerically.

Uses the negative expected Hessian of the log-likelihood.

Parameters:

• log_likelihood (callable) – Log-likelihood function L(theta).

• theta (array_like) – Parameter vector at which to evaluate.

• h (float, optional) – Step size for numerical differentiation. Default 1e-5.

Returns:

fisher – Fisher information matrix of shape (n_params, n_params).

Return type:

ndarray

Examples

>>> def log_lik(theta):
...     return -0.5 * np.sum((data - theta[0])**2 / theta[1])
>>> theta = np.array([0.0, 1.0])
>>> F = fisher_information_numerical(log_lik, theta)

Notes

For a single observation, the Fisher information is:

I(theta) = -E[d^2 log p(x|theta) / d theta^2]

This function approximates the Hessian using central differences.

pytcl.static_estimation.maximum_likelihood.fisher_information_gaussian(jacobian, noise_cov)

Fisher information for Gaussian measurement model.

For the model y = h(theta) + noise, where noise ~ N(0, R), the Fisher information is J^T R^{-1} J.

Parameters:

• jacobian (array_like) – Jacobian matrix dh/dtheta of shape (m, n).

• noise_cov (array_like) – Measurement noise covariance R of shape (m, m).

Returns:

fisher – Fisher information matrix of shape (n, n).

Return type:

ndarray

Examples

>>> H = np.array([[1, 0], [0, 1], [1, 1]])  # 3 measurements, 2 params
>>> R = np.eye(3) * 0.1
>>> F = fisher_information_gaussian(H, R)

Notes

This is the standard result for linear Gaussian models and provides the information content of measurements about parameters.

pytcl.static_estimation.maximum_likelihood.fisher_information_exponential_family(sufficient_stats, theta, data, h=1e-05)

Fisher information for exponential family distributions.

For exponential family: p(x|theta) = h(x) exp(eta(theta)^T T(x) - A(theta)) The Fisher information equals the covariance of sufficient statistics.

Parameters:

• sufficient_stats (callable) – Function T(x, theta) returning sufficient statistics.

• theta (array_like) – Natural parameters.

• data (array_like) – Observed data of shape (n_samples, …).

• h (float, optional) – Step size for numerical differentiation.

Returns:

fisher – Fisher information matrix.

Return type:

ndarray

Notes

For exponential families, I(theta) = Var[T(X)] = d^2 A(theta) / d theta^2, where A(theta) is the log-partition function.

Examples

>>> def suff_stats(x, theta):
...     return np.array([x, x**2])  # Mean and second moment
>>> data = np.random.normal(0, 1, 100)
>>> theta = np.array([0.0, 1.0])
>>> F = fisher_information_exponential_family(suff_stats, theta, data)

pytcl.static_estimation.maximum_likelihood.observed_fisher_information(log_likelihood, theta, h=1e-05)

Compute observed Fisher information (negative Hessian).

Unlike the expected Fisher information, this uses the actual observed data rather than the expectation.

Parameters:

• log_likelihood (callable) – Log-likelihood function.

• theta (array_like) – Parameter estimate.

• h (float, optional) – Step size for numerical differentiation.

Returns:

observed_fisher – Observed Fisher information matrix.

Return type:

ndarray

Notes

The observed Fisher information is often more accurate for finite samples and is asymptotically equivalent to the expected Fisher information.

Examples

>>> data = np.array([1.2, 0.8, 1.1, 0.9, 1.0])
>>> def log_lik(theta):
...     return -0.5 * np.sum((data - theta[0])**2 / theta[1]) - 2.5 * np.log(theta[1])
>>> theta = np.array([1.0, 0.1])
>>> F_obs = observed_fisher_information(log_lik, theta)

pytcl.static_estimation.maximum_likelihood.cramer_rao_bound(fisher_info, estimator_variance=None)

Compute Cramer-Rao lower bound on estimator variance.

Parameters:

• fisher_info (array_like) – Fisher information matrix of shape (n, n).

• estimator_variance (array_like, optional) – Actual estimator variance for efficiency check.

Returns:

result – Named tuple with CRB matrix, variances, and efficiency info.

Return type:

CRBResult

Examples

>>> F = np.array([[10, 0], [0, 5]])  # Fisher info
>>> result = cramer_rao_bound(F)
>>> result.variances  # Minimum achievable variances
array([0.1, 0.2])

Notes

The Cramer-Rao bound states that for any unbiased estimator:

Var(theta_hat) >= I(theta)^{-1}

The bound is achieved by efficient estimators, such as the MLE for exponential family distributions.

pytcl.static_estimation.maximum_likelihood.cramer_rao_bound_biased(fisher_info, bias_gradient)

Cramer-Rao bound for biased estimators.

Parameters:

• fisher_info (array_like) – Fisher information matrix.

• bias_gradient (array_like) – Gradient of the bias with respect to theta.

Returns:

crb – Cramer-Rao bound matrix for the biased estimator.

Return type:

ndarray

Notes

For a biased estimator with bias b(theta), the CRB becomes:

Var(theta_hat) >= (I + db/dtheta) I^{-1} (I + db/dtheta)^T

Examples

>>> F = np.array([[10.0, 0], [0, 5.0]])  # Fisher info
>>> db = np.array([[0.1, 0], [0, 0.2]])  # Bias gradient
>>> crb_biased = cramer_rao_bound_biased(F, db)
>>> crb_biased.shape
(2, 2)

pytcl.static_estimation.maximum_likelihood.efficiency(estimator_variance, crb)

Compute estimator efficiency relative to CRB.

Parameters:

• estimator_variance (array_like) – Variance of the estimator (scalar or diagonal).

• crb (array_like) – Cramer-Rao bound (scalar or diagonal).

Returns:

eff – Efficiency values in [0, 1]. Value of 1 means efficient.

Return type:

ndarray

Examples

>>> var_est = np.array([0.12, 0.25])
>>> crb = np.array([0.1, 0.2])
>>> efficiency(var_est, crb)
array([0.833, 0.8])

pytcl.static_estimation.maximum_likelihood.mle_newton_raphson(log_likelihood, score, theta_init, hessian=None, max_iter=100, tol=1e-08, h=1e-05)

Maximum likelihood estimation using Newton-Raphson.

Parameters:

• log_likelihood (callable) – Log-likelihood function L(theta).

• score (callable) – Score function (gradient of log-likelihood).

• theta_init (array_like) – Initial parameter guess.

• hessian (callable, optional) – Hessian function. If None, computed numerically.

• max_iter (int, optional) – Maximum iterations. Default 100.

• tol (float, optional) – Convergence tolerance. Default 1e-8.

• h (float, optional) – Step size for numerical Hessian.

Returns:

result – MLE result with estimate, Fisher info, and covariance.

Return type:

MLResult

Examples

>>> def log_lik(theta):
...     return -0.5 * np.sum((data - theta[0])**2)
>>> def score(theta):
...     return np.array([np.sum(data - theta[0])])
>>> result = mle_newton_raphson(log_lik, score, np.array([0.0]))

Notes

Newton-Raphson update: theta_{n+1} = theta_n - H^{-1} @ score where H is the Hessian of the log-likelihood.

pytcl.static_estimation.maximum_likelihood.mle_scoring(log_likelihood, score, fisher_info, theta_init, max_iter=100, tol=1e-08)

Maximum likelihood estimation using Fisher scoring.

Uses expected Fisher information instead of observed Hessian, which can be more stable.

Parameters:

• log_likelihood (callable) – Log-likelihood function.

• score (callable) – Score function (gradient).

• fisher_info (callable) – Fisher information function.

• theta_init (array_like) – Initial parameter guess.

• max_iter (int, optional) – Maximum iterations.

• tol (float, optional) – Convergence tolerance.

Returns:

result – MLE result.

Return type:

MLResult

Notes

Fisher scoring update: theta_{n+1} = theta_n + I(theta_n)^{-1} @ score This is equivalent to Newton-Raphson when I(theta) = -E[H].

Examples

>>> data = np.array([1.0, 1.1, 0.9, 1.2, 0.8])
>>> def log_lik(theta):
...     return -0.5 * len(data) * np.log(2*np.pi) - np.sum((data - theta[0])**2) / 2
>>> def score(theta):
...     return np.array([np.sum(data - theta[0])])
>>> def fisher(theta):
...     return np.array([[len(data)]])
>>> result = mle_scoring(log_lik, score, fisher, np.array([0.0]))

pytcl.static_estimation.maximum_likelihood.mle_gaussian(data, estimate_mean=True, estimate_variance=True)

Closed-form MLE for Gaussian distribution.

Parameters:

• data (array_like) – Observed data of shape (n_samples,) or (n_samples, n_features).

• estimate_mean (bool, optional) – Whether to estimate mean. Default True.

• estimate_variance (bool, optional) – Whether to estimate variance. Default True.

Returns:

result – MLE result with mean and/or variance estimates.

Return type:

MLResult

Examples

>>> data = np.random.normal(5, 2, 1000)
>>> result = mle_gaussian(data)
>>> result.theta  # [mean, variance]

pytcl.static_estimation.maximum_likelihood.aic(log_likelihood, n_params)

Akaike Information Criterion.

Parameters:

• log_likelihood (float) – Log-likelihood at the MLE.

• n_params (int) – Number of parameters.

Returns:

aic – AIC value (lower is better).

Return type:

float

Notes

AIC = -2 * log_likelihood + 2 * n_params

Examples

>>> log_lik = -100.0
>>> n_params = 3
>>> aic(log_lik, n_params)
206.0

pytcl.static_estimation.maximum_likelihood.bic(log_likelihood, n_params, n_samples)

Bayesian Information Criterion.

Parameters:

• log_likelihood (float) – Log-likelihood at the MLE.

• n_params (int) – Number of parameters.

• n_samples (int) – Number of samples.

Returns:

bic – BIC value (lower is better).

Return type:

float

Notes

BIC = -2 * log_likelihood + n_params * log(n_samples)

Examples

>>> log_lik = -100.0
>>> bic(log_lik, n_params=3, n_samples=100)
213.81551055796427

pytcl.static_estimation.maximum_likelihood.aicc(log_likelihood, n_params, n_samples)

Corrected Akaike Information Criterion.

Parameters:

• log_likelihood (float) – Log-likelihood at the MLE.

• n_params (int) – Number of parameters.

• n_samples (int) – Number of samples.

Returns:

aicc – AICc value (lower is better).

Return type:

float

Notes

AICc adds a correction for small sample sizes: AICc = AIC + 2*k*(k+1)/(n-k-1)

Examples

>>> log_lik = -50.0
>>> aicc(log_lik, n_params=3, n_samples=20)
109.5

Robust Estimation

M-estimators (Huber, Tukey) and RANSAC algorithms.

Robust estimation methods.

This module provides robust estimators that are resistant to outliers, including M-estimators (Huber, Tukey bisquare) and RANSAC.

References

[1]

16. 10. Huber, “Robust Statistics,” Wiley, 1981.

[2]

M. A. Fischler and R. C. Bolles, “Random Sample Consensus: A Paradigm for Model Fitting with Applications to Image Analysis and Automated Cartography,” Communications of the ACM, 1981.

class pytcl.static_estimation.robust.RobustResult(x, residuals, weights, scale, n_iter, converged)

Bases: NamedTuple

Result of robust estimation.

x

Estimated parameters.

Type:

ndarray

residuals

Residual vector.

Type:

ndarray

weights

Final weights for each observation.

Type:

ndarray

scale

Estimated scale of residuals.

Type:

float

n_iter

Number of iterations performed.

Type:

int

converged

Whether the algorithm converged.

Type:

bool

x: ndarray[tuple[Any, ...], dtype[floating]]

Alias for field number 0

residuals: ndarray[tuple[Any, ...], dtype[floating]]

Alias for field number 1

weights: ndarray[tuple[Any, ...], dtype[floating]]

Alias for field number 2

scale: float

Alias for field number 3

n_iter: int

Alias for field number 4

converged: bool

Alias for field number 5

class pytcl.static_estimation.robust.RANSACResult(x, inliers, n_inliers, residuals, n_iter, best_score)

Bases: NamedTuple

Result of RANSAC estimation.

x

Estimated parameters from best model.

Type:

ndarray

inliers

Boolean mask of inlier points.

Type:

ndarray

n_inliers

Number of inliers.

Type:

int

residuals

Residuals for all points.

Type:

ndarray

n_iter

Number of iterations performed.

Type:

int

best_score

Score of the best model found.

Type:

float

x: ndarray[tuple[Any, ...], dtype[floating]]

Alias for field number 0

inliers: ndarray[tuple[Any, ...], dtype[bool]]

Alias for field number 1

n_inliers: int

Alias for field number 2

residuals: ndarray[tuple[Any, ...], dtype[floating]]

Alias for field number 3

n_iter: int

Alias for field number 4

best_score: float

Alias for field number 5

pytcl.static_estimation.robust.huber_weight(r, c=1.345)

Huber weight function.

Parameters:

• r (array_like) – Standardized residuals.

• c (float, optional) – Tuning constant. Default 1.345 gives 95% efficiency for normal distribution.

Returns:

weights – Weights for each residual.

Return type:

ndarray

Notes

The Huber weight function is:

w(r) = 1 if |r| <= c w(r) = c / |r| if |r| > c

Examples

>>> r = np.array([0.5, 1.0, 2.0, 5.0])  # Standardized residuals
>>> w = huber_weight(r, c=1.345)
>>> w[0]  # Small residual gets weight 1
1.0
>>> w[3] < 0.5  # Large residual gets reduced weight
True

pytcl.static_estimation.robust.huber_rho(r, c=1.345)

Huber rho (loss) function.

Parameters:

• r (array_like) – Standardized residuals.

• c (float, optional) – Tuning constant.

Returns:

rho – Loss values.

Return type:

ndarray

Notes

The Huber rho function is:

rho(r) = r^2 / 2 if |r| <= c rho(r) = c * |r| - c^2/2 if |r| > c

Examples

>>> r = np.array([0.5, 1.0, 2.0])
>>> rho = huber_rho(r, c=1.345)
>>> rho[0]  # Small residual: r^2/2
0.125

pytcl.static_estimation.robust.tukey_weight(r, c=4.685)

Tukey bisquare weight function.

Parameters:

• r (array_like) – Standardized residuals.

• c (float, optional) – Tuning constant. Default 4.685 gives 95% efficiency for normal distribution.

Returns:

weights – Weights for each residual.

Return type:

ndarray

Notes

The Tukey bisquare weight function is:

w(r) = (1 - (r/c)^2)^2 if |r| <= c w(r) = 0 if |r| > c

This provides complete rejection of large outliers.

Examples

>>> r = np.array([0.5, 2.0, 5.0, 10.0])
>>> w = tukey_weight(r, c=4.685)
>>> w[0] > 0.9  # Small residual gets high weight
True
>>> w[3]  # Large residual completely rejected
0.0

pytcl.static_estimation.robust.tukey_rho(r, c=4.685)

Tukey bisquare rho (loss) function.

Parameters:

• r (array_like) – Standardized residuals.

• c (float, optional) – Tuning constant.

Returns:

rho – Loss values.

Return type:

ndarray

Notes

The Tukey rho function is:

rho(r) = c^2/6 * (1 - (1 - (r/c)^2)^3) if |r| <= c rho(r) = c^2/6 if |r| > c

Examples

>>> r = np.array([0.0, 2.0, 10.0])
>>> rho = tukey_rho(r, c=4.685)
>>> rho[0]  # Zero residual
0.0
>>> rho[2] == rho[2]  # Large residuals saturate at c^2/6
True

pytcl.static_estimation.robust.cauchy_weight(r, c=2.385)

Cauchy weight function.

Parameters:

• r (array_like) – Standardized residuals.

• c (float, optional) – Tuning constant.

Returns:

weights – Weights for each residual.

Return type:

ndarray

Notes

The Cauchy weight function is:

w(r) = 1 / (1 + (r/c)^2)

Examples

>>> r = np.array([0.0, 1.0, 5.0])
>>> w = cauchy_weight(r, c=2.385)
>>> w[0]  # Zero residual gets weight 1
1.0
>>> 0 < w[2] < 1  # Large residuals get reduced weight (but never zero)
True

pytcl.static_estimation.robust.mad(residuals, c=1.4826)

Median Absolute Deviation (MAD) scale estimator.

Parameters:

• residuals (array_like) – Residual vector.

• c (float, optional) – Consistency constant. Default 1.4826 makes MAD consistent for normal distribution.

Returns:

scale – Estimated scale.

Return type:

float

Notes

MAD = c * median(|r - median(r)|)

This is a robust scale estimator with 50% breakdown point.

Examples

>>> residuals = np.array([1.0, 1.1, 0.9, 1.0, 100.0])  # One outlier
>>> scale = mad(residuals)
>>> scale < 1.0  # Robust to the outlier
True

pytcl.static_estimation.robust.tau_scale(residuals, c1=4.5, c2=3.0)

Tau scale estimator.

Parameters:

• residuals (array_like) – Residual vector.

• c1 (float) – Tuning constants.

• c2 (float) – Tuning constants.

Returns:

scale – Estimated scale.

Return type:

float

Notes

Tau scale combines high breakdown point with efficiency.

Examples

>>> residuals = np.array([1.0, 1.1, 0.9, 1.0, 1.2, 100.0])  # One outlier
>>> scale = tau_scale(residuals)
>>> scale < 10.0  # Robust to the outlier
True

pytcl.static_estimation.robust.irls(A, b, weight_func=<function huber_weight>, scale_func=<function mad>, max_iter=50, tol=1e-06)

Iteratively Reweighted Least Squares (IRLS).

General M-estimator using IRLS algorithm.

Parameters:

• A (array_like) – Design matrix of shape (m, n).

• b (array_like) – Observation vector of shape (m,).

• weight_func (callable, optional) – Weight function w(r) for standardized residuals. Default is Huber weight.

• scale_func (callable, optional) – Scale estimation function. Default is MAD.

• max_iter (int, optional) – Maximum number of iterations. Default 50.

• tol (float, optional) – Convergence tolerance for parameter change. Default 1e-6.

Returns:

result – Robust estimation result.

Return type:

RobustResult

Examples

>>> A = np.array([[1, 1], [1, 2], [1, 3], [1, 4], [1, 10]])
>>> b = np.array([2, 3, 4, 5, 100])  # Last point is outlier
>>> result = irls(A, b)
>>> result.x  # Should be close to [1, 1] ignoring outlier

Notes

IRLS iteratively solves weighted least squares problems, updating weights based on the current residuals.

pytcl.static_estimation.robust.huber_regression(A, b, c=1.345, max_iter=50, tol=1e-06)

Huber robust regression.

M-estimation using Huber’s weight function.

Parameters:

• A (array_like) – Design matrix of shape (m, n).

• b (array_like) – Observation vector of shape (m,).

• c (float, optional) – Tuning constant. Default 1.345 gives 95% efficiency.

• max_iter (int, optional) – Maximum iterations. Default 50.

• tol (float, optional) – Convergence tolerance. Default 1e-6.

Returns:

result – Robust estimation result.

Return type:

RobustResult

Examples

>>> A = np.array([[1, 1], [1, 2], [1, 3]])
>>> b = np.array([2.1, 2.9, 100])  # Outlier in last observation
>>> result = huber_regression(A, b)

Notes

Huber regression provides a balance between efficiency for Gaussian errors and resistance to outliers.

pytcl.static_estimation.robust.tukey_regression(A, b, c=4.685, max_iter=50, tol=1e-06)

Tukey bisquare robust regression.

M-estimation using Tukey’s bisquare weight function.

Parameters:

• A (array_like) – Design matrix of shape (m, n).

• b (array_like) – Observation vector of shape (m,).

• c (float, optional) – Tuning constant. Default 4.685 gives 95% efficiency.

• max_iter (int, optional) – Maximum iterations. Default 50.

• tol (float, optional) – Convergence tolerance. Default 1e-6.

Returns:

result – Robust estimation result.

Return type:

RobustResult

Examples

>>> A = np.array([[1, 1], [1, 2], [1, 3]])
>>> b = np.array([2.1, 2.9, 100])  # Outlier in last observation
>>> result = tukey_regression(A, b)

Notes

Tukey bisquare provides complete rejection of large outliers (zero weight for residuals > c), making it more robust than Huber for gross outliers.

pytcl.static_estimation.robust.ransac(A, b, min_samples=None, residual_threshold=None, max_trials=100, stop_n_inliers=None, stop_score=None, random_state=None)

RANdom SAmple Consensus (RANSAC) regression.

Robust regression by iteratively selecting random subsets, fitting models, and identifying inliers.

Parameters:

• A (array_like) – Design matrix of shape (m, n).

• b (array_like) – Observation vector of shape (m,).

• min_samples (int, optional) – Minimum number of samples for fitting. Default is n (number of features).

• residual_threshold (float, optional) – Maximum residual for a point to be considered inlier. Default is MAD of initial OLS residuals.

• max_trials (int, optional) – Maximum number of random trials. Default 100.

• stop_n_inliers (int, optional) – Stop early if this many inliers found. Default None.

• stop_score (float, optional) – Stop early if score exceeds this. Default None.

• random_state (int, optional) – Random seed for reproducibility.

Returns:

result – RANSAC result with best model and inlier mask.

Return type:

RANSACResult

Examples

>>> rng = np.random.default_rng(42)
>>> A = np.column_stack([np.ones(100), np.arange(100)])
>>> b = 2 + 3 * np.arange(100) + rng.normal(0, 1, 100)
>>> # Add outliers
>>> b[90:] = 1000
>>> result = ransac(A, b, random_state=42)
>>> result.x  # Should be close to [2, 3]

Notes

RANSAC is particularly effective when the fraction of outliers is large (>25%), where M-estimators may struggle.

The algorithm: 1. Randomly select min_samples points 2. Fit model to selected points 3. Count inliers (points with residual < threshold) 4. If best model found, save it 5. Repeat for max_trials 6. Refit model using all inliers of best model

References

[1]

M. A. Fischler and R. C. Bolles, “Random Sample Consensus,” Communications of the ACM, 1981.

pytcl.static_estimation.robust.ransac_n_trials(n_samples, n_outliers, min_samples, probability=0.99)

Compute number of RANSAC trials needed.

Parameters:

• n_samples (int) – Total number of samples.

• n_outliers (int) – Expected number of outliers.

• min_samples (int) – Number of samples per trial.

• probability (float, optional) – Desired probability of success. Default 0.99.

Returns:

n_trials – Number of trials needed.

Return type:

int

Examples

>>> ransac_n_trials(100, 30, 2)  # 30% outliers, 2 samples per trial
10

Notes

Formula: k = log(1 - p) / log(1 - (1 - e)^n) where: - p = probability of success - e = outlier ratio - n = min_samples