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.

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

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

• 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.

Variables:

• x (ndarray) – Estimated parameters.

• residuals (ndarray) – Residual vector (y - A @ x).

• rank (int) – Effective rank of the design matrix.

• singular_values (ndarray) – Singular values of the design matrix.

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.

Variables:

• x (ndarray) – Estimated parameters.

• residuals (ndarray) – Residual vector.

• covariance (ndarray) – Estimated covariance of parameters.

• weighted_residual_sum (float) – Sum of weighted squared residuals.

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.

Variables:

• x (ndarray) – Estimated parameters.

• residuals_A (ndarray) – Corrections to the design matrix A.

• residuals_b (ndarray) – Corrections to the observation vector b.

• rank (int) – Effective rank used in the solution.

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

• 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

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

• 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.

Variables:

• theta (ndarray) – Estimated parameters.

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

• fisher_info (ndarray) – Fisher information matrix at the estimate.

• covariance (ndarray) – Estimated covariance of the parameters. This is the inverse Fisher information wherever the parameterization is unconstrained. It is not for the multivariate Gaussian, whose vec(Sigma) block carries d^2 entries for a symmetric matrix with only d(d+1)/2 free ones – there the sampling covariance is rank-deficient by construction and the inverse Fisher matrix would describe a different, unconstrained estimator (gh-20).

• n_iter (int) – Number of iterations (for iterative methods).

• converged (bool) – Whether the optimization converged.

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.

Variables:

• crb_matrix (ndarray) – Cramer-Rao bound matrix (inverse Fisher information).

• variances (ndarray) – Diagonal elements (variance bounds for each parameter).

• std_bounds (ndarray) – Square root of variances (standard deviation bounds).

• fisher_info (ndarray) – Fisher information matrix used.

• is_efficient (ndarray or None) – Boolean mask indicating which estimators achieve the bound.

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

>>> data = np.array([-1.0, 0.0, 1.0])
>>> 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)
>>> bool(np.allclose(F, [[3, 0], [0, 2]], atol=1e-4))
True

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)

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, …).

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.83333333, 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

>>> data = np.array([1.0, 2.0, 3.0])
>>> 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]))
>>> float(np.round(result.theta[0], 6))  # MLE is the sample mean
2.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

>>> rng = np.random.default_rng(42)
>>> data = rng.normal(5, 2, 1000)
>>> result = mle_gaussian(data)
>>> np.round(result.theta, 2)  # [mean, variance]
array([4.94, 3.91])
>>> bool(np.allclose(result.theta, [np.mean(data), np.var(data)]))
True

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)
107.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

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

• 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.

Variables:

• x (ndarray) – Estimated parameters.

• residuals (ndarray) – Residual vector.

• weights (ndarray) – Final weights for each observation.

• scale (float) – Estimated scale of residuals.

• n_iter (int) – Number of iterations performed.

• converged (bool) – Whether the algorithm converged.

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.

Variables:

• x (ndarray) – Estimated parameters from best model.

• inliers (ndarray) – Boolean mask of inlier points.

• n_inliers (int) – Number of inliers.

• residuals (ndarray) – Residuals for all points.

• n_iter (int) – Number of iterations performed.

• best_score (float) – Score of the best model found.

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

The MAD estimator is:

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, 5]])
>>> b = np.array([2.0, 3.0, 4.0, 100.0, 6.0])  # Fourth point is outlier
>>> result = irls(A, b)
>>> bool(np.allclose(result.x, [1, 1], atol=0.05))  # Outlier ignored
True

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)
>>> bool(np.allclose(result.x, [2, 3], atol=0.1))  # Close to [2, 3]
True

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

• 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
7

Notes

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

Localization

Closed-form static localization: TDOA least squares, bistatic range-only, range-rate-only velocity, and an ad-hoc radar covariance.

Closed-form static localization estimators.

Ports of the polynomial-free estimators from the MATLAB TCL Static_Estimation directory: TDOA least-squares emitter localization, bistatic range-only localization, range-rate-only velocity estimation, and an ad-hoc Cartesian covariance from radar sensor parameters.

References

[1]

M. D. Gillette and H. F. Silverman, “A linear closed-form algorithm for source localization from time-differences of arrival,” IEEE Signal Processing Letters, vol. 15, pp. 1-4, 2008.

[2] (1,2)

M. Malanowski and K. Kulpa, “Two methods for target localization in multistatic passive radar,” IEEE Transactions on Aerospace and Electronic Systems, vol. 48, no. 1, pp. 572-580, Jan. 2012.

[3]

D. F. Crouse, “Basic tracking using nonlinear 3D monostatic and bistatic measurements,” IEEE Aerospace and Electronic Systems Magazine, vol. 29, no. 8, Part II, pp. 4-53, Aug. 2014.

class pytcl.static_estimation.localization.DirectionOnlyLocEst(t, exit_code)

Bases: NamedTuple

Result of direction_only_static_loc_est().

Variables:

• t (ndarray) – (num_dim,) estimated target location.

• exit_code (int) – 0 on success; for algorithms 1 and 3, nonzero echoes a non-convergence status from the quasi-Newton refinement.

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

Alias for field number 0

exit_code: int

Alias for field number 1

class pytcl.static_estimation.localization.PolyStaticEst(z_cart, exit_code)

Bases: NamedTuple

Result of a polynomial-solver-based static estimator.

Variables:

• z_cart (ndarray) – (dim, num_sol) real Cartesian solutions that survived the complex and sign filters. Geometric ambiguity generally leaves more than one column.

• exit_code (int) – Exit code of poly_roots_multi_dim() (0 on success).

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

Alias for field number 0

exit_code: int

Alias for field number 1

class pytcl.static_estimation.localization.RangeOnlyLocEst(x_est, p_taylor, p_crlb)

Bases: NamedTuple

Result of range_only_static_loc_est_np().

Variables:

• x_est (ndarray) – (3,) Cartesian location estimate, or (3, 2) holding both solutions when only the minimal three measurements are given.

• p_taylor (ndarray or None) – (3, 3, num_sol) Taylor-series covariance(s), present when a measurement covariance was supplied.

• p_crlb (ndarray or None) – (3, 3, num_sol) Cramer-Rao lower bound covariance(s), present when a measurement covariance was supplied.

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

Alias for field number 0

p_taylor: ndarray[tuple[Any, ...], dtype[floating]] | None

Alias for field number 1

p_crlb: ndarray[tuple[Any, ...], dtype[floating]] | None

Alias for field number 2

pytcl.static_estimation.localization.ad_hoc_cart_cov(bandwidth, beamwidth, snr, x=None, dim=None)

Ad-hoc Cartesian covariance from radar sensor parameters.

Builds a 2D or 3D covariance whose principal axes are the range and cross-range resolutions at the estimated target location, rotated from the x-axis into the target direction.

Parameters:

• bandwidth (float) – Radar bandwidth in Hz.

• beamwidth (array_like) – Scalar beamwidth (azimuth and elevation equal), or a length-2 vector [azimuth, elevation], in radians.

• snr (float) – Signal-to-noise ratio. As in the original, the value enters the range-resolution formula directly (the MATLAB documentation calls it dB but the code applies no conversion).

• x (array_like, optional) – (2,) or (3,) estimated Cartesian target location. Default [1, 0, 0].

• dim (int, optional) – 2 for polar (range, azimuth) or 3 for spherical measurements. Default: the dimensionality of x.

Returns:

V – (dim, dim) covariance matrix.

Return type:

ndarray

Notes

Port of getAdHocCartCov.m.

Examples

>>> import numpy as np
>>> V = ad_hoc_cart_cov(5e6, [np.deg2rad(2), np.deg2rad(10)], 10.0,
...                     [1e3, 1e3, 1e3])
>>> V.shape
(3, 3)
>>> bool(np.allclose(V, V.T)) and bool(np.all(np.linalg.eigvalsh(V) > 0))
True

pytcl.static_estimation.localization.direction_only_static_loc_est(u, l_rx, algorithm=0, w=None, use_const_alg=False, r=None, r_inv=None, num_iter=1, t_init=None)

Target location from simultaneous direction (bearings) estimates.

Given unit direction vectors from at least two sensors toward a target, estimate the target’s Cartesian location in 2D or 3D.

Parameters:

• u (array_like) – (num_dim, num_meas) unit direction vectors, in the global frame, from each sensor to the target; num_meas >= 2.

• l_rx (array_like) – (num_dim, num_meas) sensor locations.

• algorithm (int, optional) –

  • 0 (default): suboptimal least-squares triangulation followed by num_iter iterations of the explicit known-range solution.

  • 1: the triangulation followed by quasi-Newton maximization of the likelihood.

  • 2: the explicit solution for known ranges r.

  • 3: quasi-Newton maximization from t_init.

• w (array_like, optional) – (num_dim, num_dim, num_meas) weights for the suboptimal triangulation (algorithms 0, 1). Default: identity.

• use_const_alg (bool, optional) – Enforce nonnegative ranges in the triangulation via a constrained solve (algorithms 0, 1). Default False.

• r (array_like, optional) – (num_meas,) known target-to-sensor ranges (algorithm 2).

• r_inv (array_like, optional) – (num_dim, num_dim, num_meas) inverse measurement covariances for the refinement stages. Default: identity. (The MATLAB original defaults these to all-ones matrices for algorithms 0 and 1 — an apparent slip; identity is used here uniformly.)

• num_iter (int, optional) – Refinement iterations for algorithms 0 and 2. Default 1.

• t_init (array_like, optional) – (num_dim,) initial estimate (algorithm 3).

Returns:

result – The location estimate and an exit code.

Return type:

DirectionOnlyLocEst

Examples

Three bearings-only sensors around a 2D target; the noise-free directions recover it.

>>> import numpy as np
>>> t_true = np.array([500.0, 800.0])
>>> l_rx = np.array([[0.0, 1000.0, -200.0], [0.0, 100.0, 900.0]])
>>> u = t_true[:, None] - l_rx
>>> u = u / np.linalg.norm(u, axis=0)
>>> res = direction_only_static_loc_est(u, l_rx)
>>> np.round(res.t, 6)
array([500., 800.])

Notes

Port of directionOnlyStaticLocEst.m, implementing the algorithms of D. F. Crouse, “Bearings-only localization using direction cosines,” Proc. 19th International Conference on Information Fusion, Jul. 2016. Documented deviations from the original, which contains several outright defects on these paths:

• triangulateKnownR honors r_inv (the original overwrites it with eye(3), discarding the weighting and crashing in 2D).

• Algorithm 2’s refinement recomputes the ranges from the current estimate t (the original recomputes them from r itself, a typo that makes its iterations meaningless).

• The quasi-Newton stages use SciPy’s BFGS with the original’s analytic gradient instead of a port of quasiNetwonBFGS (same optimum; different line-search internals).

• The MATLAB params structs are flattened into keyword arguments.

pytcl.static_estimation.localization.poly_meas_fim(x, sigma2_list, f_tx, meas_types, sensor_idx_lists, sensor_states, c=299792458.0, xi=None, w=None)

Fisher information matrix for polynomial-type localization systems.

Given simultaneous TDOA, bistatic range, emitter range-rate and/or received-frequency measurements corrupted by independent Gaussian noise, compute the Fisher information matrix (the inverse CRLB) for the location of a stationary target, using cubature integration for the expectation.

Parameters:

• x (array_like) – (num_dim,) true target location.

• sigma2_list (array_like) – (num_meas,) positive variance of each measurement.

• f_tx (float or None) – True (un-shifted) emitter frequency, required when any measurement has type 3; pass None otherwise.

• meas_types (array_like) – (num_meas,) type of each measurement: 0 TDOA, 1 bistatic range, 2 emitter range rate, 3 received frequency.

• sensor_idx_lists (array_like) – (2, num_meas) zero-based indices into sensor_states selecting the sensors of each measurement. For TDOA, row 0 is the reference sensor; for bistatic range the order does not matter; types 2 and 3 use only row 0 (set the other entry to -1 or 0).

• sensor_states (array_like) – (num_dim, num_sensors) sensor positions, or (2*num_dim, num_sensors) stacked positions and velocities (velocities are required by types 2 and 3).

• c (float, optional) – Signal propagation speed. Default: speed of light.

• xi (array_like, optional) – Cubature points (num_points, num_meas) and weights for a unit Gaussian, in pytcl’s row convention. Default: the fifth-order points for num_meas dimensions.

• w (array_like, optional) – Cubature points (num_points, num_meas) and weights for a unit Gaussian, in pytcl’s row convention. Default: the fifth-order points for num_meas dimensions.

Returns:

fim – (num_dim, num_dim) Fisher information matrix, or (num_dim+1, num_dim+1) when frequency measurements are present (the last row/column concerns the estimate of f_tx).

Return type:

ndarray

Examples

Two TDOA pairs and one bistatic range around a 3D target: the FIM is symmetric positive definite, so the position is observable.

>>> import numpy as np
>>> x = np.array([1e3, 2e3, 3e3])
>>> sensors = np.array([[0.0, 8e3, -6e3, 2e3],
...                     [0.0, 1e3, 5e3, -7e3],
...                     [0.0, -2e3, 1e3, 4e3]])
>>> idx = np.array([[0, 0, 2], [1, 2, 3]])
>>> fim = poly_meas_fim(x, [1e-14, 1e-14, 100.0], None,
...                     [0, 0, 1], idx, sensors)
>>> bool(np.all(np.linalg.eigvalsh(fim) > 0))
True

Notes

Port of computePolyMeasFIM.m, implementing the FIM equations of D. F. Crouse, “General multivariate polynomial target localization and initial estimation,” Journal of Advances in Information Fusion, vol. 13, no. 1, pp. 68-91, Jun. 2018. Unlike the MATLAB original, sensor indices are zero-based and the cubature points use pytcl’s (num_points, n) row convention.

pytcl.static_estimation.localization.range_only_static_loc_est_np(r_bi, z_loc1, z_loc2, method=1, r_cov=None)

Target location in 3D from bistatic range-only measurements.

One receiver and multiple transmitters (or vice versa); the sensors cannot all be coplanar. With noisy measurements, results degrade as the geometry approaches coplanarity.

Parameters:

• r_bi (array_like) – (num_meas,) bistatic range measurements, num_meas >= 3.

• z_loc1 (array_like) – (3, num_meas) transmitter locations (with one receiver), or receiver locations (with one transmitter).

• z_loc2 (array_like) – (3,) location of the single receiver (or transmitter). It may not be collocated with any sensor in z_loc1.

• method (int, optional) – 0 for the spherical-interpolation method of [2] (requires num_meas > 3), 1 (default) for the spherical-intersection technique of [2].

• r_cov (array_like, optional) – (num_meas, num_meas) measurement covariance. When given, the Taylor-series and CRLB covariances are computed (method 1 only, as in the original).

Returns:

result – Location estimate and, when r_cov was supplied, the two covariance estimates.

Return type:

RangeOnlyLocEst

Examples

>>> import numpy as np
>>> t = np.array([4e3, -2e3, 3e3])
>>> rx = np.array([100.0, 200.0, -50.0])
>>> tx = np.array([[0.0, 8e3, -6e3, 2e3, -3e3],
...                [0.0, 1e3, 5e3, -7e3, 2e3],
...                [0.0, -2e3, 1e3, 4e3, 9e3]])
>>> r = np.linalg.norm(t[:, None] - tx, axis=0) + np.linalg.norm(t - rx)
>>> np.round(range_only_static_loc_est_np(r, tx, rx).x_est, 6)
array([ 4000., -2000.,  3000.])

Notes

Port of rangeOnlyStaticLocEstNP.m. Two behaviors of the original are preserved deliberately: covariance outputs are only defined for method 1 (the original references variables that method 0 never creates), and the covariance of a uniquely-selected solution is linearized about solution 1’s position even when solution 2 was the one selected (the original uses xEst1 in Delta regardless of which solution won the residual comparison).

pytcl.static_estimation.localization.range_rate_ratio_to_static_pos_2d(f_rat, s_r_ref, s_rx, c=299792458.0, abs_tol=1e-09, rel_tol=1e-07, max_deg_increases=None, use_motzkin_null=False)

2D emitter location from Doppler frequency ratios alone.

A stationary emitter of UNKNOWN transmission frequency can be localized from moving sensors using only the ratios of the frequencies each sensor measures: the unknown frequency cancels in the ratio. Three sensors are needed — one reference and two others.

Parameters:

• f_rat (array_like) – (2,) frequency ratios; the numerator is the reference sensor’s measured frequency, the denominator the i-th other sensor’s.

• s_r_ref (array_like) – (4,) reference sensor state [x, y, xdot, ydot].

• s_rx (array_like) – (4, 2) states of the two other sensors.

• c (float, optional) – Propagation speed of the signal. Default: speed of light.

• abs_tol (float, optional) – Tolerances for the numerically-real and sign-consistency filters. Defaults 1e-9 and 1e-7.

• rel_tol (float, optional) – Tolerances for the numerically-real and sign-consistency filters. Defaults 1e-9 and 1e-7.

• max_deg_increases (optional) – Passed through to the polynomial solver.

• use_motzkin_null (optional) – Passed through to the polynomial solver.

Returns:

result – Real 2D solutions (the emitter plus geometric ambiguities) and the solver exit code.

Return type:

PolyStaticEst

Examples

>>> import numpy as np
>>> u_true = np.array([1e3, 5e3])
>>> ref = np.array([1000.0, 3000.0, 150.0, -150.0])
>>> s = np.array([[500.0, 1100.0], [2500.0, 2500.0]])
>>> s_dot = np.array([[300.0, 300.0], [0.0, 0.0]])
>>> c = 299792458.0
>>> rrate = lambda p, v: -v @ (u_true - p) / np.linalg.norm(u_true - p)
>>> rr_ref = rrate(ref[:2], ref[2:])
>>> f_rat = np.array(
...     [(1 - rr_ref / c) / (1 - rrate(s[:, k], s_dot[:, k]) / c) for k in (0, 1)]
... )
>>> res = range_rate_ratio_to_static_pos_2d(f_rat, ref, np.vstack([s, s_dot]))
>>> bool(np.min(np.linalg.norm(res.z_cart - u_true[:, None], axis=0)) < 1e-3)
True

Notes

Port of rangeRateRatio2StaticPos2D.m (same reference as range_rate_to_static_pos()). The problem is lifted to five variables [tx, ty, r1, r2, r3] — the target position plus the range to each sensor — before solving.

pytcl.static_estimation.localization.range_rate_to_static_pos(rr, s_rx, abs_tol=1e-09, rel_tol=1e-07, max_deg_increases=None, use_motzkin_null=False)

Stationary-emitter location from minimal range-rate measurements.

Given the minimum number of range rates for observability (2 in 2D, 3 in 3D) from moving receivers, locate a stationary emitter — e.g. drones taking Doppler measurements of a stationary phone with a known broadcast frequency. None of the receivers may be stationary.

Parameters:

• rr (array_like) – (dim,) range-rate measurements.

• s_rx (array_like) – (2*dim, dim) stacked receiver position and velocity per measurement; s_rx[:, i] = [x, y(, z), xdot, ydot(, zdot)].

• abs_tol (float, optional) – Tolerances for the numerically-real and sign-consistency filters. Defaults 1e-9 and 1e-7.

• rel_tol (float, optional) – Tolerances for the numerically-real and sign-consistency filters. Defaults 1e-9 and 1e-7.

• max_deg_increases (optional) – Passed through to the polynomial solver.

• use_motzkin_null (optional) – Passed through to the polynomial solver.

Returns:

result – Real solutions (the true emitter plus geometric ambiguities) and the solver exit code.

Return type:

PolyStaticEst

Examples

>>> import numpy as np
>>> u_true = np.array([1e3, 5e3])
>>> s = np.array([[500.0, 1100.0], [2500.0, 2500.0]])
>>> s_dot = np.array([[300.0, 300.0], [0.0, 0.0]])
>>> rr = np.array(
...     [
...         -s_dot[:, k] @ (u_true - s[:, k]) / np.linalg.norm(u_true - s[:, k])
...         for k in range(2)
...     ]
... )
>>> res = range_rate_to_static_pos(rr, np.vstack([s, s_dot]))
>>> bool(
...     np.min(np.linalg.norm(res.z_cart - u_true[:, None], axis=0)) < 1e-3
... )
True

Notes

Port of rangeRate2StaticPos.m, implementing concepts of D. F. Crouse, “General multivariate polynomial target localization and initial estimation,” Journal of Advances in Information Fusion, vol. 13, no. 1, pp. 68-91, Jun. 2018. The 3D coefficient hypermatrices are scaled by 1e-3 exactly as in the original.

pytcl.static_estimation.localization.rr_only_static_vel_est(rr, x_tx, x_rx, z_tar, use_half_range=False)

Least-squares target velocity from bistatic range-rate measurements.

Works in 2D and 3D; produces a least-squares estimate when more than the minimum number of measurements (2 in 2D, 3 in 3D) is given. Uses a non-relativistic model and ignores atmospheric effects.

Parameters:

• rr (array_like) – (num_meas,) range rates.

• x_tx (array_like or None) – (2*d, num_meas) stacked transmitter position/velocity states, or a single (2*d,) state shared by all measurements. Pass None when the target itself is the transmitter (an emitter).

• x_rx (array_like) – (2*d, num_meas) stacked receiver states, or a single (2*d,) state shared by all measurements.

• z_tar (array_like) – (d,) Cartesian target position.

• use_half_range (bool, optional) – True if the range rates are one-way (monostatic convention); they are doubled internally. Default False.

Returns:

v_est – (d,) least-squares Cartesian velocity estimate.

Return type:

ndarray

Examples

An emitter (the target is the transmitter) observed by three moving receivers; error-free one-way range rates recover its velocity:

>>> import numpy as np
>>> z_tar = np.array([1.5, -0.4, 2.2])
>>> v_tar = np.array([0.3, 1.1, -0.7])
>>> x_rx = np.array([[0.5, -1.2, 2.0],
...                  [1.0, 0.3, -1.5],
...                  [-0.6, 1.8, 0.4],
...                  [0.1, -0.5, 0.7],
...                  [-0.2, 0.4, 0.1],
...                  [0.3, 0.2, -0.4]])
>>> h = z_tar[:, None] - x_rx[:3]
>>> h = h / np.linalg.norm(h, axis=0)
>>> rr = np.sum(h * (v_tar[:, None] - x_rx[3:]), axis=0)
>>> np.round(rr_only_static_vel_est(rr, None, x_rx, z_tar), 9)
array([ 0.3,  1.1, -0.7])

Notes

Port of RROnlyStaticVelEst.m, implementing Equation 41 in Section IV E of [3], with the target-is-transmitter case handled specially to remove the singularity.

pytcl.static_estimation.localization.tdoa_only_static_loc_est(time_delays, ref_rx_locs, non_ref_rx_locs, c=299792458.0)

Closed-form least-squares emitter location from TDOA measurements.

A minimum of one reference receiver and four TDOA measurements is needed for observability in 3D. For minimal (exactly-determined) systems use tdoa_to_cart instead (not yet ported).

Parameters:

• time_delays (array_like or sequence of array_like) – With a single reference receiver, an (n,) vector of time differences between each receiver and the reference. With multiple references, a sequence whose i-th element holds the delay vector for the receivers paired with the i-th reference. The form (array or sequence of arrays) must match non_ref_rx_locs.

• ref_rx_locs (array_like) – (3,) location of the single reference receiver, or (3, num_refs) locations of all reference receivers.

• non_ref_rx_locs (array_like or sequence of array_like) – With a single reference, a (3, n) matrix of receiver locations. With multiple references, a sequence whose i-th element is the (3, n_i) matrix of receivers paired with the i-th reference.

• c (float, optional) – Propagation speed. Default: speed of light.

Returns:

source_loc – (3,) emitter location. Exact in an error-free setting; otherwise a least-squares solution with respect to a non-standard cost function.

Return type:

ndarray

Examples

>>> import numpy as np
>>> t = np.array([27.0, 0.0, -42.0])
>>> ref = np.array([9.0, 39.0, 100.0])
>>> rx = np.array([[65.0, 64.0, -128.0, 0.0],
...                [10.0, 71.0, 6.0, -20.0],
...                [-60.0, 43.0, 12.0, 4.0]])
>>> c = 341.0
>>> tdoa = (np.linalg.norm(t[:, None] - rx, axis=0)
...         - np.linalg.norm(t - ref)) / c
>>> np.round(tdoa_only_static_loc_est(tdoa, ref, rx, c), 9) + 0.0
array([ 27.,   0., -42.])

Notes

Port of TDOAOnlyStaticLocEst.m, implementing the linear closed-form algorithm of [1].

pytcl.static_estimation.localization.tdoa_to_cart(tdoa, l_rx1, l_rx2, c=299792458.0, abs_tol=1e-09, rel_tol=1e-07, max_deg_increases=None, use_motzkin_null=False)

Target location from a minimal set of TDOA measurements.

Exactly 2 measurements in 2D or 3 in 3D — the minimal number for observability, unlike the overdetermined least-squares tdoa_only_static_loc_est(). The hyperbolic equations are turned into simultaneous multivariate polynomials and solved with poly_roots_multi_dim().

Parameters:

• tdoa (array_like) – (dim,) time differences of arrival; tdoa[i] is the arrival time at l_rx2[:, i] minus the arrival time at l_rx1[:, i].

• l_rx1 (array_like) – (dim, dim) reference sensor positions, one column per measurement, or a single (dim,) position shared by all measurements.

• l_rx2 (array_like) – (dim, dim) non-reference sensor positions.

• c (float, optional) – Propagation speed. Default: speed of light.

• abs_tol (float, optional) – Tolerances used both to decide whether a root is numerically real and to discard sign-flipped ghost solutions introduced by the squaring. Defaults 1e-9 and 1e-7.

• rel_tol (float, optional) – Tolerances used both to decide whether a root is numerically real and to discard sign-flipped ghost solutions introduced by the squaring. Defaults 1e-9 and 1e-7.

• max_deg_increases (int, optional) – Passed through to the polynomial solver.

• use_motzkin_null (bool, optional) – Passed through to the polynomial solver. Default False.

Returns:

result – Real solutions and the solver exit code. As in the original, if the sign filter would discard every candidate, the first one is kept anyway.

Return type:

PolyStaticEst

Examples

>>> import numpy as np
>>> S1 = np.array([9.0, 39.0, 100.0])
>>> S2 = np.array([65.0, 10.0, -60.0])
>>> S3 = np.array([64.0, 71.0, 43.0])
>>> S4 = np.array([-128.0, 6.0, 12.0])
>>> t = np.array([27.0, 0.0, -42.0])
>>> c = 341.0
>>> d = lambda a, b: np.linalg.norm(t - a) - np.linalg.norm(t - b)
>>> tdoa = np.array([d(S2, S1), d(S3, S1), d(S3, S4)]) / c
>>> l_rx1 = np.column_stack([S1, S1, S4])
>>> l_rx2 = np.column_stack([S2, S3, S3])
>>> res = tdoa_to_cart(tdoa, l_rx1, l_rx2, c)
>>> np.round(res.z_cart[:, 0], 6) + 0.0  # +0.0 normalizes signed zeros
array([ 27.,   0., -42.])

Notes

Port of TDOA2Cart.m, following the polynomial formulation of M. P. Williams, “Solving polynomial equations using linear algebra,” Johns Hopkins Technical Digest, vol. 28, no. 4, pp. 354-363, 2010 (with the sign-of-u typo of the paper fixed, as in the original).