Static Estimation

This example demonstrates static estimation algorithms for parameter estimation, sensor calibration, and model fitting.

Overview

Static estimation methods are fundamental for:

• Sensor calibration: Bias and scale factor estimation

• Model fitting: Parameter estimation from data

• Regression analysis: Relationship between variables

• Data fusion: Combining multiple measurements

Least Squares Methods

Ordinary Least Squares (OLS)

• Minimizes sum of squared residuals

• Assumes equal measurement uncertainty

• Optimal for Gaussian noise

Weighted Least Squares (WLS)

• Accounts for varying measurement precision

• Weights = 1/variance for optimal results

• Essential for heteroscedastic data

Total Least Squares (TLS)

• Errors in both dependent and independent variables

• Also known as errors-in-variables regression

• Avoids bias from noisy predictors

Generalized Least Squares (GLS)

• Accounts for correlated measurement errors

• Uses full covariance matrix

Recursive Least Squares (RLS)

• Online/sequential estimation

• Updates estimate with each new measurement

• Useful for real-time applications

Ridge Regression

• L2 regularization for ill-conditioned problems

• Shrinks coefficients toward zero

• Handles multicollinearity

Robust Estimation

Huber M-estimator

• Combines L2 (small residuals) and L1 (large residuals)

• Less sensitive to outliers than OLS

• Iteratively reweighted least squares (IRLS)

Tukey Bisquare M-estimator

• Completely rejects large outliers (zero weight)

• More aggressive outlier handling

• Identifies outliers through weights

RANSAC

• Random Sample Consensus

• Robust to high outlier percentages

• Separates inliers from outliers

Maximum Likelihood

MLE for Gaussian

• Estimates mean and variance

• Optimal for Gaussian data

Fisher Information

• Measures information in data about parameters

• Determines estimation precision limits

Cramer-Rao Bound

• Lower bound on estimator variance

• Efficiency = CRB / actual variance

Model Selection

AIC (Akaike Information Criterion)

• Balances fit quality and model complexity

• Penalizes additional parameters

• Lower is better

BIC (Bayesian Information Criterion)

• Stronger penalty for complexity

• Consistent model selection

• Prefers simpler models than AIC

Code Highlights

The example demonstrates:

• OLS with ordinary_least_squares()

• WLS with weighted_least_squares()

• TLS with total_least_squares()

• RLS with recursive_least_squares()

• Ridge with ridge_regression()

• Huber with huber_regression()

• Tukey with tukey_regression()

• RANSAC with ransac()

• MLE with mle_gaussian()

• Model selection with aic(), bic()

Source Code

"""
Static Estimation Example
=========================

This example demonstrates static estimation algorithms in PyTCL:

Least Squares Methods:
- Ordinary Least Squares (OLS)
- Weighted Least Squares (WLS)
- Total Least Squares (TLS)
- Generalized Least Squares (GLS)
- Recursive Least Squares (RLS)
- Ridge Regression

Robust Estimation:
- Huber M-estimator
- Tukey bisquare M-estimator
- RANSAC for outlier-robust fitting

Maximum Likelihood Estimation:
- MLE for Gaussian parameters
- Fisher Information and Cramer-Rao Bounds
- Model selection (AIC, BIC)

These methods are fundamental for parameter estimation, sensor calibration,
and model fitting in the presence of noise and outliers.
"""

from pathlib import Path

import numpy as np
import plotly.graph_objects as go
from plotly.subplots import make_subplots

# Output directory for generated plots
OUTPUT_DIR = Path(__file__).parent.parent / "docs" / "_static" / "images" / "examples"
OUTPUT_DIR.mkdir(parents=True, exist_ok=True)

# Global flag to control plotting
SHOW_PLOTS = True


from pytcl.static_estimation import (  # Least squares; Robust estimation; MLE and Fisher Information; Model selection
    aic,
    aicc,
    bic,
    cramer_rao_bound,
    efficiency,
    fisher_information_gaussian,
    fisher_information_numerical,
    generalized_least_squares,
    huber_regression,
    irls,
    mad,
    mle_gaussian,
    ordinary_least_squares,
    ransac,
    recursive_least_squares,
    ridge_regression,
    total_least_squares,
    tukey_regression,
    weighted_least_squares,
)


def demo_ordinary_least_squares():
    """Demonstrate ordinary least squares."""
    print("=" * 70)
    print("Ordinary Least Squares Demo")
    print("=" * 70)

    np.random.seed(42)

    # Generate linear regression data
    n_samples = 50
    x = np.linspace(0, 10, n_samples)
    true_slope = 2.5
    true_intercept = 1.0
    noise_std = 0.5

    y = true_intercept + true_slope * x + np.random.randn(n_samples) * noise_std

    # Design matrix [1, x]
    A = np.column_stack([np.ones(n_samples), x])

    # OLS solution
    result = ordinary_least_squares(A, y)

    print(f"\nTrue parameters: intercept={true_intercept}, slope={true_slope}")
    print(f"OLS estimate: intercept={result.x[0]:.4f}, slope={result.x[1]:.4f}")
    print(f"\nResidual sum of squares: {np.sum(result.residuals**2):.4f}")
    print(f"Matrix rank: {result.rank}")

    # Coefficient of determination
    ss_res = np.sum(result.residuals**2)
    ss_tot = np.sum((y - np.mean(y)) ** 2)
    r_squared = 1 - ss_res / ss_tot
    print(f"R² = {r_squared:.4f}")

    # Plot OLS fit
    if SHOW_PLOTS:
        fig = make_subplots(
            rows=1,
            cols=2,
            subplot_titles=[
                f"Ordinary Least Squares (R² = {r_squared:.4f})",
                "Residual Plot",
            ],
        )

        # Fit plot
        fig.add_trace(
            go.Scatter(
                x=x,
                y=y,
                mode="markers",
                marker=dict(color="blue", size=8, opacity=0.6),
                name="Data",
            ),
            row=1,
            col=1,
        )

        x_line = np.linspace(x.min(), x.max(), 100)
        y_true_line = true_intercept + true_slope * x_line
        y_fit_line = result.x[0] + result.x[1] * x_line

        fig.add_trace(
            go.Scatter(
                x=x_line,
                y=y_true_line,
                mode="lines",
                line=dict(color="green", width=2, dash="dash"),
                name="True line",
            ),
            row=1,
            col=1,
        )

        fig.add_trace(
            go.Scatter(
                x=x_line,
                y=y_fit_line,
                mode="lines",
                line=dict(color="red", width=2),
                name="OLS fit",
            ),
            row=1,
            col=1,
        )

        # Residuals plot
        residuals = y - (result.x[0] + result.x[1] * x)
        fig.add_trace(
            go.Scatter(
                x=x,
                y=residuals,
                mode="markers",
                marker=dict(color="blue", size=8, opacity=0.6),
                name="Residuals",
                showlegend=False,
            ),
            row=1,
            col=2,
        )

        fig.add_hline(y=0, line_dash="dash", line_color="red", row=1, col=2)

        fig.update_xaxes(title_text="x", row=1, col=1)
        fig.update_yaxes(title_text="y", row=1, col=1)
        fig.update_xaxes(title_text="x", row=1, col=2)
        fig.update_yaxes(title_text="Residual", row=1, col=2)

        fig.update_layout(height=500, width=1000, showlegend=True)
        fig.write_html(str(OUTPUT_DIR / "static_ols.html"))
        print("\n  [Plot saved to static_ols.html]")


def demo_weighted_least_squares():
    """Demonstrate weighted least squares."""
    print("\n" + "=" * 70)
    print("Weighted Least Squares Demo")
    print("=" * 70)

    np.random.seed(42)

    # Data with heteroscedastic noise (varying variance)
    n_samples = 50
    x = np.linspace(0, 10, n_samples)
    true_slope = 2.0
    true_intercept = 3.0

    # Noise increases with x
    noise_std = 0.2 + 0.1 * x
    y = true_intercept + true_slope * x + np.random.randn(n_samples) * noise_std

    A = np.column_stack([np.ones(n_samples), x])

    # OLS (ignores varying noise)
    result_ols = ordinary_least_squares(A, y)

    # WLS with weights = 1/variance
    weights = 1 / noise_std**2
    result_wls = weighted_least_squares(A, y, weights=weights)

    print(f"\nTrue parameters: intercept={true_intercept}, slope={true_slope}")
    print(
        f"\nOLS estimate: intercept={result_ols.x[0]:.4f}, slope={result_ols.x[1]:.4f}"
    )
    print(f"WLS estimate: intercept={result_wls.x[0]:.4f}, slope={result_wls.x[1]:.4f}")

    print("\nNote: WLS gives more weight to precise measurements (low variance)")
    print("and typically produces better estimates when noise is heteroscedastic.")


def demo_total_least_squares():
    """Demonstrate total least squares (errors-in-variables)."""
    print("\n" + "=" * 70)
    print("Total Least Squares Demo")
    print("=" * 70)

    np.random.seed(42)

    # Both x and y have measurement errors
    n_samples = 30
    x_true = np.linspace(0, 10, n_samples)
    true_slope = 1.5
    true_intercept = 2.0

    # Add noise to both x and y
    x_noise_std = 0.3
    y_noise_std = 0.5

    x = x_true + np.random.randn(n_samples) * x_noise_std
    y = true_intercept + true_slope * x_true + np.random.randn(n_samples) * y_noise_std

    A = np.column_stack([np.ones(n_samples), x])

    # OLS (assumes x is error-free)
    result_ols = ordinary_least_squares(A, y)

    # TLS (accounts for errors in x)
    result_tls = total_least_squares(A, y)

    print(f"\nTrue parameters: intercept={true_intercept}, slope={true_slope}")
    print(
        f"\nOLS estimate: intercept={result_ols.x[0]:.4f}, slope={result_ols.x[1]:.4f}"
    )
    print(f"TLS estimate: intercept={result_tls.x[0]:.4f}, slope={result_tls.x[1]:.4f}")

    print("\nNote: TLS is preferred when independent variables have measurement error.")
    print("OLS typically underestimates the true slope in this case.")


def demo_recursive_least_squares():
    """Demonstrate recursive least squares for online estimation."""
    print("\n" + "=" * 70)
    print("Recursive Least Squares Demo")
    print("=" * 70)

    np.random.seed(42)

    # Online parameter estimation
    n_samples = 100
    x = np.linspace(0, 10, n_samples)
    true_slope = 2.0
    true_intercept = 1.0
    noise_std = 0.3

    y = true_intercept + true_slope * x + np.random.randn(n_samples) * noise_std

    # Initialize RLS (2 parameters: intercept and slope)
    n_params = 2
    x_est = np.zeros(n_params)  # Initial parameter estimate
    P = np.eye(n_params) * 100.0  # Initial covariance (large uncertainty)

    print("\nOnline parameter estimation with RLS:")
    print("-" * 50)

    checkpoints = [10, 25, 50, 100]
    checkpoint_idx = 0

    for i in range(n_samples):
        # Measurement vector [1, x_i] for y = intercept + slope * x
        a = np.array([1.0, x[i]])
        y_i = y[i]

        # RLS update
        x_est, P = recursive_least_squares(x_est, P, a, y_i)

        # Print at checkpoints
        if checkpoint_idx < len(checkpoints) and i + 1 == checkpoints[checkpoint_idx]:
            print(
                f"  After {i+1:>3} samples: intercept={x_est[0]:.4f}, "
                f"slope={x_est[1]:.4f}"
            )
            checkpoint_idx += 1

    print(f"\nTrue values: intercept={true_intercept}, slope={true_slope}")
    print("\nNote: RLS converges to true values as more data arrives.")


def demo_ridge_regression():
    """Demonstrate ridge regression for ill-conditioned problems."""
    print("\n" + "=" * 70)
    print("Ridge Regression Demo")
    print("=" * 70)

    np.random.seed(42)

    # Create collinear predictors (ill-conditioned)
    n_samples = 30
    x1 = np.random.randn(n_samples)
    x2 = x1 + np.random.randn(n_samples) * 0.1  # x2 ≈ x1

    true_coef = np.array([1.0, 2.0, 3.0])  # [intercept, b1, b2]
    y = true_coef[0] + true_coef[1] * x1 + true_coef[2] * x2
    y += np.random.randn(n_samples) * 0.5

    A = np.column_stack([np.ones(n_samples), x1, x2])

    # OLS solution (may be unstable)
    result_ols = ordinary_least_squares(A, y)

    # Ridge regression with regularization
    lambdas = [0.0, 0.01, 0.1, 1.0]

    print(f"\nTrue coefficients: {true_coef}")
    print("\nEstimates with different regularization:")
    print("-" * 60)
    print(f"{'Lambda':>10} {'Intercept':>12} {'b1':>12} {'b2':>12}")
    print("-" * 60)

    for lam in lambdas:
        if lam == 0:
            x_hat = result_ols.x
        else:
            x_hat = ridge_regression(A, y, alpha=lam)  # Returns array directly
        print(f"{lam:>10.2f} {x_hat[0]:>12.4f} {x_hat[1]:>12.4f} " f"{x_hat[2]:>12.4f}")

    print("\nNote: Ridge regression shrinks coefficients toward zero,")
    print("which helps stabilize estimates for collinear predictors.")


def demo_robust_estimation():
    """Demonstrate robust estimation methods."""
    print("\n" + "=" * 70)
    print("Robust Estimation Demo")
    print("=" * 70)

    np.random.seed(42)

    # Generate data with outliers
    n_samples = 50
    n_outliers = 5

    x = np.linspace(0, 10, n_samples)
    true_slope = 2.0
    true_intercept = 1.0

    y = true_intercept + true_slope * x + np.random.randn(n_samples) * 0.5

    # Add outliers
    outlier_idx = np.random.choice(n_samples, n_outliers, replace=False)
    y[outlier_idx] += np.random.randn(n_outliers) * 10

    A = np.column_stack([np.ones(n_samples), x])

    print(f"\nData: {n_samples} samples with {n_outliers} outliers")
    print(f"True parameters: intercept={true_intercept}, slope={true_slope}")

    # OLS (sensitive to outliers)
    result_ols = ordinary_least_squares(A, y)
    print(f"\nOLS: intercept={result_ols.x[0]:.4f}, slope={result_ols.x[1]:.4f}")

    # Huber M-estimator
    result_huber = huber_regression(A, y)
    print(f"Huber: intercept={result_huber.x[0]:.4f}, slope={result_huber.x[1]:.4f}")
    print(f"  Iterations: {result_huber.n_iter}, Converged: {result_huber.converged}")

    # Tukey bisquare M-estimator
    result_tukey = tukey_regression(A, y)
    print(f"Tukey: intercept={result_tukey.x[0]:.4f}, slope={result_tukey.x[1]:.4f}")
    print(f"  Iterations: {result_tukey.n_iter}, Converged: {result_tukey.converged}")

    # Analyze weights from Tukey estimator
    weights = result_tukey.weights
    detected_outliers = np.where(weights < 0.1)[0]
    print(f"\nDetected outliers (low weight): {detected_outliers}")
    print(f"True outlier indices: {sorted(outlier_idx)}")

    # Plot robust estimation comparison
    if SHOW_PLOTS:
        fig = make_subplots(
            rows=1,
            cols=2,
            subplot_titles=[
                "Robust Estimation: OLS vs M-estimators",
                "Tukey M-estimator Weights (red = true outliers)",
            ],
        )

        # Fit comparison
        # Regular data points
        regular_mask = np.ones(n_samples, dtype=bool)
        regular_mask[outlier_idx] = False

        fig.add_trace(
            go.Scatter(
                x=x[regular_mask],
                y=y[regular_mask],
                mode="markers",
                marker=dict(color="blue", size=8, opacity=0.6),
                name="Data",
            ),
            row=1,
            col=1,
        )

        fig.add_trace(
            go.Scatter(
                x=x[outlier_idx],
                y=y[outlier_idx],
                mode="markers",
                marker=dict(color="red", size=10, opacity=0.8),
                name="Outliers",
            ),
            row=1,
            col=1,
        )

        x_line = np.linspace(x.min(), x.max(), 100)

        fig.add_trace(
            go.Scatter(
                x=x_line,
                y=true_intercept + true_slope * x_line,
                mode="lines",
                line=dict(color="green", width=2, dash="dash"),
                name="True",
            ),
            row=1,
            col=1,
        )

        fig.add_trace(
            go.Scatter(
                x=x_line,
                y=result_ols.x[0] + result_ols.x[1] * x_line,
                mode="lines",
                line=dict(color="black", width=2),
                name="OLS",
            ),
            row=1,
            col=1,
        )

        fig.add_trace(
            go.Scatter(
                x=x_line,
                y=result_huber.x[0] + result_huber.x[1] * x_line,
                mode="lines",
                line=dict(color="magenta", width=2, dash="dash"),
                name="Huber",
            ),
            row=1,
            col=1,
        )

        fig.add_trace(
            go.Scatter(
                x=x_line,
                y=result_tukey.x[0] + result_tukey.x[1] * x_line,
                mode="lines",
                line=dict(color="cyan", width=2, dash="dot"),
                name="Tukey",
            ),
            row=1,
            col=1,
        )

        # Weights from Tukey estimator
        colors = ["red" if i in outlier_idx else "blue" for i in range(n_samples)]
        fig.add_trace(
            go.Scatter(
                x=x,
                y=weights,
                mode="markers",
                marker=dict(color=colors, size=8, opacity=0.7),
                name="Weights",
                showlegend=False,
            ),
            row=1,
            col=2,
        )

        fig.add_hline(y=0.1, line_dash="dash", line_color="red", row=1, col=2)

        fig.update_xaxes(title_text="x", row=1, col=1)
        fig.update_yaxes(title_text="y", row=1, col=1)
        fig.update_xaxes(title_text="x", row=1, col=2)
        fig.update_yaxes(title_text="Tukey weight", row=1, col=2)

        fig.update_layout(height=500, width=1200, showlegend=True)
        fig.write_html(str(OUTPUT_DIR / "static_robust_estimation.html"))
        print("\n  [Plot saved to static_robust_estimation.html]")


def demo_ransac():
    """Demonstrate RANSAC for robust fitting."""
    print("\n" + "=" * 70)
    print("RANSAC Demo")
    print("=" * 70)

    np.random.seed(42)

    # Line fitting with many outliers
    n_inliers = 40
    n_outliers = 20

    # Inliers: points on line y = 2x + 1
    x_inliers = np.random.uniform(0, 10, n_inliers)
    y_inliers = 1 + 2 * x_inliers + np.random.randn(n_inliers) * 0.3

    # Outliers: random points
    x_outliers = np.random.uniform(0, 10, n_outliers)
    y_outliers = np.random.uniform(-5, 25, n_outliers)

    x = np.concatenate([x_inliers, x_outliers])
    y = np.concatenate([y_inliers, y_outliers])

    # Shuffle
    perm = np.random.permutation(len(x))
    x, y = x[perm], y[perm]

    A = np.column_stack([np.ones(len(x)), x])

    print(f"\nData: {n_inliers} inliers + {n_outliers} outliers = {len(x)} total")
    print("True line: y = 2x + 1")

    # OLS
    result_ols = ordinary_least_squares(A, y)
    print(f"\nOLS: y = {result_ols.x[1]:.4f}x + {result_ols.x[0]:.4f}")

    # RANSAC
    threshold = 1.0  # Inlier threshold (residual threshold)

    result_ransac = ransac(
        A, y, min_samples=2, residual_threshold=threshold, max_trials=100
    )

    print(f"RANSAC: y = {result_ransac.x[1]:.4f}x + {result_ransac.x[0]:.4f}")
    print(f"  Inliers found: {result_ransac.n_inliers}/{len(x)}")

    print("\nNote: RANSAC successfully recovers the true line despite")
    print("33% of the data being outliers.")

    # Plot RANSAC
    if SHOW_PLOTS:
        fig = go.Figure()

        # Determine inliers based on RANSAC residuals
        y_pred = result_ransac.x[0] + result_ransac.x[1] * x
        residuals = np.abs(y - y_pred)
        is_inlier = residuals < threshold

        fig.add_trace(
            go.Scatter(
                x=x[is_inlier],
                y=y[is_inlier],
                mode="markers",
                marker=dict(color="blue", size=8, opacity=0.6),
                name="Inliers",
            )
        )

        fig.add_trace(
            go.Scatter(
                x=x[~is_inlier],
                y=y[~is_inlier],
                mode="markers",
                marker=dict(color="red", size=8, opacity=0.6),
                name="Outliers",
            )
        )

        x_line = np.linspace(x.min(), x.max(), 100)

        fig.add_trace(
            go.Scatter(
                x=x_line,
                y=1 + 2 * x_line,
                mode="lines",
                line=dict(color="green", width=2, dash="dash"),
                name="True line",
            )
        )

        fig.add_trace(
            go.Scatter(
                x=x_line,
                y=result_ols.x[0] + result_ols.x[1] * x_line,
                mode="lines",
                line=dict(color="black", width=2),
                name="OLS",
            )
        )

        fig.add_trace(
            go.Scatter(
                x=x_line,
                y=result_ransac.x[0] + result_ransac.x[1] * x_line,
                mode="lines",
                line=dict(color="red", width=2),
                name="RANSAC",
            )
        )

        fig.update_layout(
            title=f"RANSAC Line Fitting ({result_ransac.n_inliers} inliers / {len(x)} total)",
            xaxis_title="x",
            yaxis_title="y",
            height=500,
            width=800,
            showlegend=True,
        )
        fig.write_html(str(OUTPUT_DIR / "static_ransac.html"))
        print("\n  [Plot saved to static_ransac.html]")


def demo_mle_and_fisher():
    """Demonstrate MLE and Fisher information."""
    print("\n" + "=" * 70)
    print("Maximum Likelihood Estimation Demo")
    print("=" * 70)

    np.random.seed(42)

    # Generate Gaussian data
    n_samples = 100
    true_mean = 5.0
    true_std = 2.0

    data = true_mean + np.random.randn(n_samples) * true_std

    print(f"\nTrue parameters: mean={true_mean}, std={true_std}")
    print(f"Sample size: {n_samples}")

    # MLE for Gaussian parameters
    # theta[0] = mean, theta[1] = variance
    result = mle_gaussian(data)
    mle_mean = result.theta[0]
    mle_var = result.theta[1]
    mle_std = np.sqrt(mle_var)
    print(f"\nMLE estimates: mean={mle_mean:.4f}, std={mle_std:.4f}")
    print(f"Log-likelihood: {result.log_likelihood:.4f}")

    # Fisher information
    # For Gaussian: I(μ) = n/σ², I(σ²) = n/(2σ⁴)
    fisher_mean = n_samples / true_std**2
    fisher_var = n_samples / (2 * true_std**4)

    print(f"\nFisher information:")
    print(f"  I(mean) = {fisher_mean:.4f}")
    print(f"  I(variance) = {fisher_var:.4f}")

    # Cramer-Rao bounds
    crb_mean = 1 / fisher_mean
    crb_var = 1 / fisher_var

    print(f"\nCramer-Rao lower bounds (variance of unbiased estimator):")
    print(f"  Var(mean_hat) ≥ {crb_mean:.6f}")
    print(f"  Var(var_hat) ≥ {crb_var:.6f}")

    # Check actual MLE variance (through simulation)
    n_trials = 1000
    mean_estimates = []
    for _ in range(n_trials):
        sample = true_mean + np.random.randn(n_samples) * true_std
        result = mle_gaussian(sample)
        mean_estimates.append(result.theta[0])

    actual_variance = np.var(mean_estimates)
    print(f"\nSimulated variance of mean estimator: {actual_variance:.6f}")
    print(f"Efficiency: {crb_mean / actual_variance:.4f}")
    print("(Efficiency = 1.0 means estimator achieves the CRB)")


def demo_model_selection():
    """Demonstrate model selection using information criteria."""
    print("\n" + "=" * 70)
    print("Model Selection Demo")
    print("=" * 70)

    np.random.seed(42)

    # True model: quadratic y = 1 + 2x + 0.5x²
    n_samples = 50
    x = np.linspace(0, 5, n_samples)
    y_true = 1 + 2 * x + 0.5 * x**2
    y = y_true + np.random.randn(n_samples) * 0.5

    print("\nTrue model: y = 1 + 2x + 0.5x²")
    print("Comparing polynomial models of degree 1 through 5:")
    print("-" * 60)
    print(f"{'Degree':>6} {'RSS':>10} {'k':>4} {'AIC':>10} {'BIC':>10}")
    print("-" * 60)

    results = []
    for degree in range(1, 6):
        # Build design matrix for polynomial
        A = np.column_stack([x**i for i in range(degree + 1)])
        result = ordinary_least_squares(A, y)

        rss = np.sum(result.residuals**2)
        k = degree + 1  # Number of parameters
        n = n_samples

        # Compute log-likelihood (assuming Gaussian errors)
        sigma2 = rss / n
        log_lik = -n / 2 * np.log(2 * np.pi * sigma2) - rss / (2 * sigma2)

        aic_val = aic(log_lik, k)
        bic_val = bic(log_lik, k, n)

        results.append((degree, rss, k, aic_val, bic_val))
        print(f"{degree:>6} {rss:>10.4f} {k:>4} {aic_val:>10.2f} {bic_val:>10.2f}")

    # Find best models
    best_aic = min(results, key=lambda x: x[3])
    best_bic = min(results, key=lambda x: x[4])

    print(f"\nBest by AIC: degree {best_aic[0]}")
    print(f"Best by BIC: degree {best_bic[0]}")
    print("\nNote: True model is degree 2. AIC/BIC help avoid overfitting.")

    # Plot model selection
    if SHOW_PLOTS:
        fig = make_subplots(
            rows=1,
            cols=2,
            subplot_titles=["Polynomial Model Fits", "Model Selection: AIC vs BIC"],
        )

        # Model fits
        fig.add_trace(
            go.Scatter(
                x=x,
                y=y,
                mode="markers",
                marker=dict(color="blue", size=8, opacity=0.6),
                name="Data",
            ),
            row=1,
            col=1,
        )

        x_line = np.linspace(x.min(), x.max(), 100)
        fig.add_trace(
            go.Scatter(
                x=x_line,
                y=1 + 2 * x_line + 0.5 * x_line**2,
                mode="lines",
                line=dict(color="green", width=2, dash="dash"),
                name="True (deg 2)",
            ),
            row=1,
            col=1,
        )

        # Fit degrees 1, 2, 3
        colors = ["red", "green", "purple"]
        for i, deg in enumerate([1, 2, 3]):
            A_fit = np.column_stack([x_line**j for j in range(deg + 1)])
            A_data = np.column_stack([x**j for j in range(deg + 1)])
            coef = ordinary_least_squares(A_data, y).x
            y_fit = A_fit @ coef
            fig.add_trace(
                go.Scatter(
                    x=x_line,
                    y=y_fit,
                    mode="lines",
                    line=dict(width=1.5, dash="solid" if deg == 2 else "dot"),
                    name=f"Degree {deg}",
                ),
                row=1,
                col=1,
            )

        # AIC/BIC comparison
        degrees = [r[0] for r in results]
        aic_vals = [r[3] for r in results]
        bic_vals = [r[4] for r in results]

        fig.add_trace(
            go.Bar(
                x=[d - 0.2 for d in degrees],
                y=aic_vals,
                width=0.35,
                name="AIC",
                marker_color="blue",
                opacity=0.7,
            ),
            row=1,
            col=2,
        )

        fig.add_trace(
            go.Bar(
                x=[d + 0.2 for d in degrees],
                y=bic_vals,
                width=0.35,
                name="BIC",
                marker_color="orange",
                opacity=0.7,
            ),
            row=1,
            col=2,
        )

        fig.update_xaxes(title_text="x", row=1, col=1)
        fig.update_yaxes(title_text="y", row=1, col=1)
        fig.update_xaxes(title_text="Polynomial Degree", row=1, col=2)
        fig.update_yaxes(title_text="Information Criterion", row=1, col=2)

        fig.update_layout(height=500, width=1200, showlegend=True)
        fig.write_html(str(OUTPUT_DIR / "static_model_selection.html"))
        print("\n  [Plot saved to static_model_selection.html]")


def main():
    """Run all demonstrations."""
    print("\n" + "#" * 70)
    print("# PyTCL Static Estimation Example")
    print("#" * 70)

    # Least squares methods
    demo_ordinary_least_squares()
    demo_weighted_least_squares()
    demo_total_least_squares()
    demo_recursive_least_squares()
    demo_ridge_regression()

    # Robust estimation
    demo_robust_estimation()
    demo_ransac()

    # MLE and model selection
    demo_mle_and_fisher()
    demo_model_selection()

    print("\n" + "=" * 70)
    print("Example complete!")
    if SHOW_PLOTS:
        print("Plots saved: static_ols.html, static_robust_estimation.html,")
        print("             static_ransac.html, static_model_selection.html")
    print("=" * 70)


if __name__ == "__main__":
    main()

Running the Example

python examples/static_estimation.py

See Also

• Kalman Filter Comparison - Dynamic estimation

• performance_evaluation - Tracking metrics

• gaussian_mixtures - Clustering and mixture models