Robust Regression with Huber Loss (IRLS)

Least squares lets a few outliers drag the whole fit. The Huber loss is quadratic for small residuals and linear for large ones:

L_\delta(r) = \begin{cases}\tfrac12 r^2 & |r| \le \delta\\ \delta(|r| - \tfrac12\delta) & \text{otherwise}\end{cases}

Fit it with iteratively reweighted least squares: start from the OLS solution, then repeat: compute residuals, set weights w = 1 where |r| <= delta and delta / |r| elsewhere, and solve the weighted least-squares problem. Stop when the coefficients change by less than tol.

Implement:

huber_loss(r, delta): the mean loss
fit_huber(X, y, delta=1.0, iters=100, tol=1e-8): with an intercept. Return (intercept, coefs).
fit_ols(X, y): the same return format, for comparison