Robust Estimation Tutorial

This tutorial demonstrates robust regression methods for handling outliers and non-Gaussian noise, comparing three estimators on the same outlier-contaminated linear dataset.

Topics covered:

  • Ordinary Least Squares (OLS) as the outlier-sensitive baseline

  • RANSAC (Random Sample Consensus)

  • Iteratively Reweighted Least Squares (IRLS) with a Huber weight function

  • Comparing residual distributions across estimators

Contaminated Data

The true model is y = 2x + 1; 15% of the 100 points are corrupted with large uniform outliers on top of the baseline Gaussian noise:

import numpy as np

np.random.seed(42)
x_data = np.linspace(0, 10, 100)
y_data = 2 * x_data + 1 + np.random.randn(100) * 0.5

outlier_idx = np.random.choice(100, 15, replace=False)
y_data[outlier_idx] += np.random.uniform(-5, 5, 15)

Ordinary Least Squares

A = np.vstack([x_data, np.ones_like(x_data)]).T
params_ols = np.linalg.lstsq(A, y_data, rcond=None)[0]

OLS minimizes squared residuals over all points, so the 15% outlier fraction pulls the fitted slope and intercept measurably off the true model.

RANSAC

RANSAC repeatedly fits a line to a minimal 2-point sample, keeps the sample whose fit has the most inliers under a residual threshold, and reports the inlier-only fit:

best_inliers = 0
for _ in range(100):
    sample_idx = np.random.choice(100, 2, replace=False)
    params = np.linalg.lstsq(
        np.vstack([x_data[sample_idx], np.ones(2)]).T, y_data[sample_idx],
        rcond=None,
    )[0]
    inlier_mask = np.abs(y_data - (params[0] * x_data + params[1])) < 1.0
    if inlier_mask.sum() > best_inliers:
        best_inliers, best_params, best_inlier_mask = inlier_mask.sum(), params, inlier_mask

Because it never lets an outlier influence the fit once it is excluded as an outlier, RANSAC recovers a slope/intercept close to the true model even with 15% contamination.

Iteratively Reweighted Least Squares

IRLS starts from the OLS fit and iteratively down-weights points with large residuals using a Huber weight (w = 1 inside the threshold, w = c / |residual| outside it), refitting a weighted least-squares problem each iteration:

params = params_ols.copy()
for _ in range(10):
    residuals = y_data - (params[0] * x_data + params[1])
    weights = np.where(np.abs(residuals) <= 1.0, 1.0, 1.0 / np.abs(residuals))
    W = np.diag(weights)
    A = np.vstack([x_data, np.ones_like(x_data)]).T
    params = np.linalg.lstsq(A.T @ W @ A, A.T @ W @ y_data, rcond=None)[0]

Unlike RANSAC, IRLS never fully discards a point – even a downweighted outlier still contributes a small residual pull – so its fit sits between OLS and RANSAC in practice.

Next Steps