(Robust) Huber Regression - Namir - 2026-03-01 18:44
I asked Claude AI to write a function huberMLR() to implement the robust Huber Regression algorithm that handles outliers very well. The function handles the matrix x (with n rows and m columns)and the column vector y *with n rows), The function returns the regression coefficients, adjusted R-suare, and the square root of sum of squared errors. I also asked Claude AI for a test program,
Here is the code for function huberMLR():
Code: // ============================================================
// huberMLR() -- Huber Robust Multiple Linear Regression
// for the HP Prime Graphing Calculator
// ============================================================
//
// Method
// ------
// Huber regression minimises the Huber loss:
//
// L(r) = 0.5 * r^2 if |r| <= delta
// L(r) = delta*(|r| - 0.5*d) if |r| > delta
//
// This makes the estimator robust to outliers: residuals
// within delta are penalised quadratically (like OLS) and
// residuals beyond delta are penalised linearly.
//
// The tuning parameter delta is estimated automatically:
// sigma = MAD / 0.6745
// delta = 1.345 * sigma
// where MAD = median of |residuals|. The factor 1.345
// gives 95% asymptotic efficiency relative to OLS on
// clean Gaussian data.
//
// Algorithm
// ---------
// Iteratively Reweighted Least Squares (IRLS):
// 1. Initialise with OLS: b = (X'X)^-1 X'y
// 2. Compute residuals r = y - X*b
// 3. Estimate sigma from MAD(|r|)
// 4. Set delta = 1.345 * sigma
// 5. Compute Huber weights:
// w_i = 1 if |r_i| <= delta
// w_i = delta / |r_i| if |r_i| > delta
// 6. Solve WLS: b = (X'WX)^-1 X'Wy
// 7. Repeat from 2 until convergence
//
// To avoid creating an n x n weight matrix, X'WX and X'Wy
// are accumulated row-by-row in O(n*(m+1)^2) space.
//
// Parameters
// ----------
// x : data matrix (n rows, m columns)
// y : response column vector (n rows, 1 column)
//
// Returns
// -------
// { coeffs, adjR2, rootSSE }
//
// coeffs : row vector [intercept, b1, b2, ..., bm]
// (1 row, m+1 columns)
//
// adjR2 : weighted adjusted R-squared using final
// Huber weights. Outlier points are down-
// weighted so that adjR2 reflects the fit
// quality on the clean portion of data.
// wSSE = sum( w_i*(y_i - yHat_i)^2 )
// wSST = sum( w_i*(y_i - yBarW)^2 )
// adjR2 = 1 - (wSSE/(n-p)) / (wSST/(n-1))
//
// rootSSE : sqrt( sum( (y_i - yHat_i)^2 ) )
// where yHat is computed from x and the
// regression coefficients.
//
// Usage
// -----
// x := [[1],[2],[3],[4],[5]];
// y := [[3],[5],[7],[9],[11]];
// huberMLR(x, y) --> { [[1,2]], adjR2, rootSSE }
//
// ============================================================
EXPORT huberMLR(x, y)
BEGIN
LOCAL n, m, p, i, j, k;
LOCAL maxIt, tol, iter;
LOCAL xA, b, bOld;
LOCAL r, rAbs, w;
LOCAL med, sig, delta;
LOCAL xWx, xWy, diff;
LOCAL yHat;
LOCAL SSE;
LOCAL adjR2, rootSSE;
LOCAL coeffs;
LOCAL L, half, lst;
LOCAL sumW, yBarW, wSSE, wSST;
lst := SIZE(x);
n := lst(1);
m := lst(2);
p := m + 1; // number of parameters
maxIt := 50;
tol := 1E-10;
// ========================================================
// BUILD AUGMENTED MATRIX xA = [1 | x]
// xA is n x p; first column = 1 (intercept)
// ========================================================
xA := MAKEMAT(0, n, p);
FOR i FROM 1 TO n DO
xA(i, 1) := 1;
FOR j FROM 1 TO m DO
xA(i, j + 1) := x(i, j);
END;
END;
// ========================================================
// INITIAL OLS ESTIMATE b = (X'X)^-1 X'y
// ========================================================
b := (TRN(xA) * xA)^(-1) * TRN(xA) * y;
// ========================================================
// IRLS LOOP
// ========================================================
FOR iter FROM 1 TO maxIt DO
bOld := b;
// ---- Residuals ----
r := y - xA * b;
// ---- Estimate sigma via MAD ----
// Build sorted list of |residuals|
L := {};
FOR i FROM 1 TO n DO
L := CONCAT(L, {ABS(r(i, 1))});
END;
L := SORT(L);
// Median of |r|
IF n MOD 2 == 1 THEN
med := L((n + 1) / 2);
ELSE
half := n / 2;
med := (L(half) + L(half + 1)) / 2;
END;
// Scale estimate and tuning parameter
sig := med / 0.6745;
IF sig < 1E-12 THEN
sig := 1E-12;
END;
delta := 1.345 * sig;
// ---- Build X'WX and X'Wy row by row ----
xWx := MAKEMAT(0, p, p);
xWy := MAKEMAT(0, p, 1);
FOR i FROM 1 TO n DO
// Huber weight
rAbs := ABS(r(i, 1));
IF rAbs <= delta THEN
w := 1;
ELSE
w := delta / rAbs;
END;
// Accumulate weighted outer product
FOR j FROM 1 TO p DO
xWy(j, 1) := xWy(j, 1)
+ w * xA(i, j) * y(i, 1);
FOR k FROM 1 TO p DO
xWx(j, k) := xWx(j, k)
+ w * xA(i, j) * xA(i, k);
END;
END;
END;
// ---- Solve WLS ----
b := xWx^(-1) * xWy;
// ---- Convergence check ----
diff := 0;
FOR i FROM 1 TO p DO
diff := diff + (b(i, 1) - bOld(i, 1))^2;
END;
diff := √(diff);
IF diff < tol THEN
BREAK;
END;
END; // iter
// ========================================================
// COMPUTE STATISTICS
// ========================================================
// Fitted values and final residuals
yHat := xA * b;
r := y - yHat;
// ---- Recompute final Huber weights from converged b ----
L := {};
FOR i FROM 1 TO n DO
L := CONCAT(L, {ABS(r(i, 1))});
END;
L := SORT(L);
IF n MOD 2 == 1 THEN
med := L((n + 1) / 2);
ELSE
half := n / 2;
med := (L(half) + L(half + 1)) / 2;
END;
sig := med / 0.6745;
IF sig < 1E-12 THEN
sig := 1E-12;
END;
delta := 1.345 * sig;
// ---- Weighted mean of y ----
sumW := 0;
yBarW := 0;
FOR i FROM 1 TO n DO
rAbs := ABS(r(i, 1));
IF rAbs <= delta THEN
w := 1;
ELSE
w := delta / rAbs;
END;
sumW := sumW + w;
yBarW := yBarW + w * y(i, 1);
END;
yBarW := yBarW / sumW;
// ---- Weighted SSE and SST ----
// wSSE = sum( w_i * (y_i - yHat_i)^2 )
// wSST = sum( w_i * (y_i - yBarW)^2 )
// Outlier points get small w, so they no longer
// dominate the sums.
wSSE := 0;
wSST := 0;
FOR i FROM 1 TO n DO
rAbs := ABS(r(i, 1));
IF rAbs <= delta THEN
w := 1;
ELSE
w := delta / rAbs;
END;
wSSE := wSSE + w * (y(i, 1) - yHat(i, 1))^2;
wSST := wSST + w * (y(i, 1) - yBarW)^2;
END;
// ---- Weighted adjusted R-squared ----
IF wSST > 0 AND n > p THEN
adjR2 := 1 - (wSSE / (n - p))
/ (wSST / (n - 1));
ELSE
adjR2 := 0;
END;
// ---- Root SSE (unweighted, on all points) ----
SSE := 0;
FOR i FROM 1 TO n DO
SSE := SSE + (y(i, 1) - yHat(i, 1))^2;
END;
rootSSE := √(SSE);
// Coefficients as a row vector
coeffs := TRN(b);
// ========================================================
// RETURN { coeffs, adjR2, rootSSE }
// ========================================================
RETURN {coeffs, adjR2, rootSSE};
END;
and here is the code for testHuber:
Code: // ============================================================
// testHuber() -- Test harness for huberMLR()
// for the HP Prime Graphing Calculator
// ============================================================
//
// Runs 5 test cases that exercise the Huber regression on
// clean data, data with outliers, and multi-predictor data.
// For cases with outliers, OLS coefficients are shown beside
// the Huber coefficients so the robustness is visible.
//
// Test Cases
// ----------
// 1. y = 2x + 3 (clean, no outliers)
// 2. y = 2x + 3 (same but 2 outliers added)
// 3. y = 3x1 + 2x2 + 1 (clean, 2 predictors)
// 4. y = 3x1 + 2x2 + 1 (same with 2 outliers)
// 5. Larger dataset with ~20% outliers
//
// Usage
// -----
// Store both huberMLR and testHuber on the calculator,
// then run testHuber().
//
// ============================================================
// ---- Helper: OLS for comparison ----------------------------
olsCoeffs(x, y)
BEGIN
LOCAL n, m, p, i, j;
LOCAL xA, b, lst;
lst := SIZE(x);
n := lst(1);
m := lst(2);;
p := m + 1;
xA := MAKEMAT(0, n, p);
FOR i FROM 1 TO n DO
xA(i, 1) := 1;
FOR j FROM 1 TO m DO
xA(i, j + 1) := x(i, j);
END;
END;
b := (TRN(xA) * xA)^(-1) * TRN(xA) * y;
RETURN TRN(b);
END;
// ---- Helper: format coefficient vector as string -----------
fmtVec(v)
BEGIN
LOCAL s, j, c, lst;
lst := SIZE(v);
c := lst(2);
s := "[";
FOR j FROM 1 TO c DO
IF j > 1 THEN
s := s + ", ";
END;
s := s + STRING(v(1, j), 2, 4);
END;
s := s + "]";
RETURN s;
END;
// ============================================================
// MAIN TEST PROGRAM
// ============================================================
EXPORT testHuber()
BEGIN
LOCAL x, y, xv, result;
LOCAL coeffs, adjR2, rootSSE;
LOCAL olsC;
LOCAL nPass, nFail;
LOCAL i, j, v;
LOCAL tol;
nPass := 0;
nFail := 0;
PRINT();
PRINT("================================");
PRINT(" huberMLR() Test Suite");
PRINT("================================");
PRINT("");
// ========================================================
// TEST 1: y = 2x + 3 (clean data, no outliers)
// True: intercept=3, slope=2
// ========================================================
PRINT("TEST 1: y=2x+3 (clean)");
PRINT(" True coeffs: [3, 2]");
x := [[1],[2],[3],[4],[5],
[6],[7],[8],[9],[10]];
y := MAKEMAT(0, 10, 1);
FOR i FROM 1 TO 10 DO
y(i, 1) := 2 * x(i, 1) + 3;
END;
result := huberMLR(x, y);
coeffs := result(1);
adjR2 := result(2);
rootSSE := result(3);
PRINT(" Huber : " + fmtVec(coeffs));
PRINT(" adjR2 : " + STRING(adjR2, 2, 6));
PRINT(" rtSSE : " + STRING(rootSSE, 2, 6));
// Check intercept ~ 3, slope ~ 2
tol := 0.1;
IF ABS(coeffs(1,1)-3)<tol
AND ABS(coeffs(1,2)-2)<tol THEN
PRINT(" --> PASS");
nPass := nPass + 1;
ELSE
PRINT(" --> FAIL");
nFail := nFail + 1;
END;
PRINT("");
// ========================================================
// TEST 2: y = 2x + 3 with 2 outliers
// Points 3 and 7 get large y-outliers.
// Huber should still recover ~[3, 2].
// OLS will be distorted.
// ========================================================
PRINT("TEST 2: y=2x+3 + 2 outliers");
PRINT(" True coeffs: [3, 2]");
y(3, 1) := 200; // outlier at x=3
y(7, 1) := -150; // outlier at x=7
result := huberMLR(x, y);
coeffs := result(1);
adjR2 := result(2);
rootSSE := result(3);
olsC := olsCoeffs(x, y);
PRINT(" OLS : " + fmtVec(olsC));
PRINT(" Huber : " + fmtVec(coeffs));
PRINT(" adjR2 : " + STRING(adjR2, 2, 6));
PRINT(" rtSSE : " + STRING(rootSSE, 2, 6));
// Huber slope should be closer to 2 than OLS
IF ABS(coeffs(1,2) - 2)
< ABS(olsC(1,2) - 2) THEN
PRINT(" Huber closer to true slope");
PRINT(" --> PASS");
nPass := nPass + 1;
ELSE
PRINT(" --> FAIL");
nFail := nFail + 1;
END;
PRINT("");
// ========================================================
// TEST 3: y = 3*x1 + 2*x2 + 1 (clean, 2 predictors)
// True: [1, 3, 2]
// ========================================================
PRINT("TEST 3: y=3x1+2x2+1 (clean)");
PRINT(" True coeffs: [1, 3, 2]");
x := MAKEMAT(0, 12, 2);
y := MAKEMAT(0, 12, 1);
v := 0;
FOR i FROM 1 TO 4 DO
FOR j FROM 1 TO 3 DO
v := v + 1;
x(v, 1) := i;
x(v, 2) := j;
y(v, 1) := 3 * i + 2 * j + 1;
END;
END;
result := huberMLR(x, y);
coeffs := result(1);
adjR2 := result(2);
rootSSE := result(3);
PRINT(" Huber : " + fmtVec(coeffs));
PRINT(" adjR2 : " + STRING(adjR2, 2, 6));
PRINT(" rtSSE : " + STRING(rootSSE, 2, 6));
tol := 0.1;
IF ABS(coeffs(1,1)-1)<tol
AND ABS(coeffs(1,2)-3)<tol
AND ABS(coeffs(1,3)-2)<tol THEN
PRINT(" --> PASS");
nPass := nPass + 1;
ELSE
PRINT(" --> FAIL");
nFail := nFail + 1;
END;
PRINT("");
// ========================================================
// TEST 4: y = 3*x1 + 2*x2 + 1 with outliers
// Corrupt 2 of 12 points.
// ========================================================
PRINT("TEST 4: y=3x1+2x2+1 + outliers");
PRINT(" True coeffs: [1, 3, 2]");
y(2, 1) := 500; // outlier
y(10, 1) := -300; // outlier
result := huberMLR(x, y);
coeffs := result(1);
adjR2 := result(2);
rootSSE := result(3);
olsC := olsCoeffs(x, y);
PRINT(" OLS : " + fmtVec(olsC));
PRINT(" Huber : " + fmtVec(coeffs));
PRINT(" adjR2 : " + STRING(adjR2, 2, 6));
PRINT(" rtSSE : " + STRING(rootSSE, 2, 6));
// Check Huber slope1 closer to 3 than OLS
IF ABS(coeffs(1,2) - 3)
< ABS(olsC(1,2) - 3) THEN
PRINT(" Huber closer to true slope");
PRINT(" --> PASS");
nPass := nPass + 1;
ELSE
PRINT(" --> FAIL");
nFail := nFail + 1;
END;
PRINT("");
// ========================================================
// TEST 5: Larger dataset, ~20% outliers
// y = 5*x + 10, n=20, 4 outliers
// ========================================================
PRINT("TEST 5: y=5x+10, n=20, 4 outliers");
PRINT(" True coeffs: [10, 5]");
x := MAKEMAT(0, 20, 1);
y := MAKEMAT(0, 20, 1);
FOR i FROM 1 TO 20 DO
x(i, 1) := i;
y(i, 1) := 5 * i + 10;
END;
// Add 4 outliers (20%)
y(4, 1) := 800;
y(9, 1) := -600;
y(14, 1) := 900;
y(18, 1) := -500;
result := huberMLR(x, y);
coeffs := result(1);
adjR2 := result(2);
rootSSE := result(3);
olsC := olsCoeffs(x, y);
PRINT(" OLS : " + fmtVec(olsC));
PRINT(" Huber : " + fmtVec(coeffs));
PRINT(" adjR2 : " + STRING(adjR2, 2, 6));
PRINT(" rtSSE : " + STRING(rootSSE, 2, 6));
// Huber slope should be closer to 5
IF ABS(coeffs(1,2) - 5)
< ABS(olsC(1,2) - 5) THEN
PRINT(" Huber closer to true slope");
PRINT(" --> PASS");
nPass := nPass + 1;
ELSE
PRINT(" --> FAIL");
nFail := nFail + 1;
END;
PRINT("");
// ========================================================
// SUMMARY
// ========================================================
PRINT("================================");
PRINT(" RESULTS: "
+ STRING(nPass) + " passed, "
+ STRING(nFail) + " failed (of 5)");
PRINT("================================");
RETURN {nPass, nFail};
END;
|