Least Squares with Outliers

Least squares gives every data point the same say in the fit, and it measures disagreement by the square of the residual. A measurement that is wrong by a large amount therefore pulls the fit toward itself with a force that grows with its error. This section measures how far one bad point can move a least squares line, and why its position along the

x

axis matters as much as its error. It then fits the same data three other ways: with robust losses solved by iteratively reweighted least squares, with a Student’s t noise model, and with the inlier and outlier mixture of the EM algorithm section. The last part fits a model of the Texas electricity load to a winter that contains Winter Storm Uri, and shows which of these methods keep the storm hours from bending the fit. The notation is that of the linear regression section: pairs

(x^{(i)}, y^{(i)})

for

i = 1, \ldots, m

, basis functions

\boldsymbol\phi(x)

, parameters

\boldsymbol\theta

, the design matrix

\boldsymbol\Phi

with rows

\boldsymbol\phi(x^{(i)})^\top

, and Gaussian noise of variance

\sigma^2

. Equation numbers refer to Bishop (Bishop, n.d.), who writes

\mathbf w

for

\boldsymbol\theta

t_n

for

y^{(i)}

N

for

m

, and a precision

\beta = 1/\sigma^2

import numpy as np
import pandas as pd
from scipy.stats import norm

The running example

The data are the line

y = 1 + 0.5x

observed at 60 points with Gaussian noise of standard deviation

\sigma = 0.3

, the same generator as in the EM section. The model is the straight line

g(x; \boldsymbol\theta) = \boldsymbol\theta^\top \boldsymbol\phi(x)

with

\boldsymbol\phi(x) = (1, x)^\top

, so

\theta_0

is the intercept and

\theta_1

the slope. Least squares minimizes the sum-of-squares error

E(\boldsymbol\theta) = \frac{1}{2}\sum_{i=1}^{m} e_i^2, \qquad e_i = y^{(i)} - \boldsymbol\theta^\top \boldsymbol\phi(x^{(i)}),

which is Bishop’s Eq. 3.12. Here

e_i

is the residual of point

i

. The minimizer solves the normal equations,

\boldsymbol\theta = (\boldsymbol\Phi^\top \boldsymbol\Phi)^{-1} \boldsymbol\Phi^\top \mathbf y

(Eq. 3.15). It is also the maximum likelihood estimate when the noise is Gaussian, as the conditional maximum likelihood section shows.

def make_line_data(m=60, n_out=0, seed=1, sigma=0.3):
    # y = 1 + 0.5 x + noise; the last n_out targets are replaced by uniform gross errors
    rng = np.random.default_rng(seed)
    x = rng.uniform(0, 10, m)
    y = 1.0 + 0.5 * x + sigma * rng.standard_normal(m)
    is_outlier = np.zeros(m, dtype=bool)
    if n_out:
        is_outlier[m - n_out:] = True
        y[is_outlier] = rng.uniform(-4, 12, n_out)
    return x, y, is_outlier


def design(x):
    return np.column_stack([np.ones_like(x), x])


def least_squares(Phi, y):
    return np.linalg.lstsq(Phi, y, rcond=None)[0]


theta_true = np.array([1.0, 0.5])
x_clean, y_clean, _ = make_line_data()
Phi_clean = design(x_clean)
print("least squares on clean data: theta =", np.round(least_squares(Phi_clean, y_clean), 3))

least squares on clean data: theta = [0.929 0.509]

One bad point and leverage

Add a 61st point to the clean data and move its target

y

away from the line by an amount

\Delta

. Each target enters the least squares solution linearly, because

\boldsymbol\theta = (\boldsymbol\Phi^\top \boldsymbol\Phi)^{-1} \boldsymbol\Phi^\top \mathbf y

is a fixed matrix times

\mathbf y

. The estimate therefore moves in proportion to

\Delta

, with no limit, so a single bad point among 61 can move the slope by any amount. In the language of robust statistics, the breakdown point of least squares, the smallest fraction of bad points that can carry the estimate arbitrarily far, is

1/m

, which goes to zero as

m

grows. How far it moves depends on where the point sits. The fitted values are

\hat{\mathbf y} = \mathbf H \mathbf y

with the hat matrix

\mathbf H = \boldsymbol\Phi(\boldsymbol\Phi^\top \boldsymbol\Phi)^{-1}\boldsymbol\Phi^\top

. Its diagonal element

h_{ii} = \boldsymbol\phi(x^{(i)})^\top(\boldsymbol\Phi^\top\boldsymbol\Phi)^{-1}\boldsymbol\phi(x^{(i)})

is the leverage of point

i

: the amount by which the fitted value at

x^{(i)}

moves when

y^{(i)}

moves by one unit. Points near the mean of the

x^{(i)}

have low leverage. Points at the ends of the range, or beyond it, have high leverage. The code places the bad point at

x = 5

, near the middle of the data, and at

x = 15

, beyond its end, and fits for a range of

\Delta

def add_point(x, y, x_new, y_new):
    return np.append(x, x_new), np.append(y, y_new)


deltas = np.linspace(0, 40, 41)
slopes_ls = {}
for x_bad in (5.0, 15.0):
    xs, ys = add_point(x_clean, y_clean, x_bad, 1.0 + 0.5 * x_bad)
    Phi_b = design(xs)
    h = np.diag(Phi_b @ np.linalg.inv(Phi_b.T @ Phi_b) @ Phi_b.T)
    slopes_ls[x_bad] = np.array([least_squares(Phi_b, ys + d * (np.arange(len(ys)) == len(ys) - 1))[1] for d in deltas])
    rate = (slopes_ls[x_bad][-1] - slopes_ls[x_bad][0]) / deltas[-1]
    print(f"bad point at x = {x_bad:4.1f}: leverage h = {h[-1]:.3f} (mean leverage {h.mean():.3f}), "
          f"slope changes by {rate:+.4f} per unit of Delta")

bad point at x =  5.0: leverage h = 0.017 (mean leverage 0.033), slope changes by -0.0006 per unit of Delta
bad point at x = 15.0: leverage h = 0.176 (mean leverage 0.033), slope changes by +0.0167 per unit of Delta

Least squares slope against the error added to one point. The line for a point at x = 15 climbs steeply and without limit; the line for a point at x = 5 stays close to the true slope of 0.5.

Least squares slope as one point of 61 is moved away from the line, for a point near the middle of the data and a point beyond its end. A point near the middle has almost no lever arm on the slope: moving it shifts the intercept, and the slope barely changes. The same error at

x = 15

tilts the line. Leverage depends only on the inputs

x^{(i)}

, not on the targets, so a bad target at a high-leverage input does the most damage, and it is also the hardest to detect from the residuals. The fit moves toward the point, so the point’s own residual stays small. The grid example at the end of this section shows the same effect.

Robust losses and iteratively reweighted least squares

A robust loss replaces the square with a function

\rho

that grows more slowly for large residuals:

E_\rho(\boldsymbol\theta) = \sum_{i=1}^{m} \rho\!\left(\frac{e_i}{s}\right),

where

s

is a scale for the residuals, close to

\sigma

for the inliers. Two common choices are the Huber loss, which is quadratic up to a threshold

\delta

and linear beyond it, and the Cauchy loss, which grows only logarithmically:

\rho_{\text{Huber}}(a) = \begin{cases} \tfrac12 a^2 & |a| \le \delta \\ \delta |a| - \tfrac12 \delta^2 & |a| > \delta \end{cases} \qquad \rho_{\text{Cauchy}}(a) = \frac{c^2}{2}\ln\!\left(1 + \frac{a^2}{c^2}\right).

Setting the gradient of

E_\rho

to zero gives

\sum_i \rho'(e_i/s)\,\boldsymbol\phi(x^{(i)}) = 0

. Write

\rho'(a) = w(a)\,a

, with the weight function

w(a) = \rho'(a)/a

. The condition becomes

\sum_i w_i\, e_i\, \boldsymbol\phi(x^{(i)}) = 0

with

w_i = w(e_i/s)

, which is the normal equation of a weighted least squares problem:

\boldsymbol\theta = (\boldsymbol\Phi^\top \mathbf W \boldsymbol\Phi)^{-1} \boldsymbol\Phi^\top \mathbf W \mathbf y, \qquad \mathbf W = \operatorname{diag}(w_1, \ldots, w_m).

This has the form of Bishop’s Eq. 14.42. The weights depend on the residuals, which depend on

\boldsymbol\theta

, so the solution is found by iteration: compute the residuals, compute the weights, solve the weighted problem, repeat. This is iteratively reweighted least squares (IRLS). The weight functions are

w_{\text{Huber}}(a) = \min\!\left(1, \frac{\delta}{|a|}\right), \qquad w_{\text{Cauchy}}(a) = \frac{1}{1 + a^2/c^2}.

Huber weights fall off as

1/|a|

, so a point’s influence on the fit,

w(a)\,a

, stops growing beyond

\delta

but never shrinks. Cauchy weights fall off as

1/a^2

, so the influence of a far point goes back toward zero. The scan registration section uses the same Huber weights inside ICP. The code estimates the scale at every iteration from the median absolute deviation of the residuals,

s = 1.4826 \cdot \operatorname{median}_i |e_i - \operatorname{median}_j e_j|

, which equals

\sigma

for Gaussian residuals and ignores up to half the points. The constants

\delta = 1.345

and

c = 2.385

are the usual ones: with Gaussian noise and no outliers they lose only 5% of the efficiency of least squares.

def mad_scale(e):
    # median absolute deviation, scaled to equal sigma for Gaussian residuals
    return 1.4826 * np.median(np.abs(e - np.median(e)))


def huber_weight(a, delta=1.345):
    return np.minimum(1.0, delta / np.maximum(np.abs(a), 1e-12))


def cauchy_weight(a, c=2.385):
    return 1.0 / (1.0 + (a / c) ** 2)


def weighted_least_squares(Phi, y, w):
    sw = np.sqrt(w)
    return np.linalg.lstsq(Phi * sw[:, None], y * sw, rcond=None)[0]


def irls(Phi, y, weight, theta0, n_iter=100):
    # iteratively reweighted least squares with the scale re-estimated each iteration
    theta = theta0
    for _ in range(n_iter):
        e = y - Phi @ theta
        w = weight(e / mad_scale(e))
        theta = weighted_least_squares(Phi, y, w)
    return theta, w

Repeat the single-point experiment with the two robust losses, for the point at

x = 15

slopes_robust = {"Huber": [], "Cauchy": []}
xs, ys0 = add_point(x_clean, y_clean, 15.0, 1.0 + 0.5 * 15.0)
Phi_b = design(xs)
for d in deltas:
    ys = ys0 + d * (np.arange(len(ys0)) == len(ys0) - 1)
    theta_ls = least_squares(Phi_b, ys)
    slopes_robust["Huber"].append(irls(Phi_b, ys, huber_weight, theta_ls)[0][1])
    slopes_robust["Cauchy"].append(irls(Phi_b, ys, cauchy_weight, theta_ls)[0][1])
for name, s in slopes_robust.items():
    print(f"{name:7s} slope at Delta = 0: {s[0]:.3f}   largest slope over all Delta: {max(s):.3f}   at Delta = 40: {s[-1]:.3f}")
print(f"least squares slope at Delta = 40: {slopes_ls[15.0][-1]:.3f}")

Huber   slope at Delta = 0: 0.500   largest slope over all Delta: 0.515   at Delta = 40: 0.515
Cauchy  slope at Delta = 0: 0.502   largest slope over all Delta: 0.514   at Delta = 40: 0.504
least squares slope at Delta = 40: 1.177

Slope against the error added to the high-leverage point for three methods. Least squares rises steadily off the top of the plot. Huber rises a little and levels off. Cauchy rises slightly for small errors and then returns to the true slope.

Slope estimates as the high-leverage point moves away from the line. Huber caps the influence of the point, so the slope stops moving once the residual passes the threshold, but it settles at a biased value. Cauchy weights the point less and less as it moves away, so the slope returns to the true value. The Cauchy loss is not convex, so the weighted problem can have several solutions, and IRLS finds the one nearest its start. The Huber loss is convex, so IRLS reaches the same answer from any start.

The Student’s t view

Robust losses can also be read as likelihoods. Bishop’s Section 2.3.7 introduces the Student’s t distribution

\mathrm{St}(x \mid \mu, \lambda, \nu)

(Eq. 2.159) and shows in Figure 2.16 that a maximum likelihood fit under it is much less affected by a few outliers than a Gaussian fit. The reason is its form as a scale mixture of Gaussians (Eq. 2.160):

\mathrm{St}(x \mid \mu, \lambda, \nu) = \int_0^\infty \mathcal N\big(x \mid \mu, (\eta\lambda)^{-1}\big)\, \mathrm{Gam}(\eta \mid \nu/2, \nu/2)\, d\eta.

Each observation has its own precision scale

\eta

, drawn from a gamma distribution, so some points are allowed to be far noisier than others. Take the t distribution as the noise model of the regression,

p(y^{(i)} \mid x^{(i)}) = \mathrm{St}(y^{(i)} \mid \boldsymbol\theta^\top\boldsymbol\phi(x^{(i)}), \lambda, \nu)

. Its negative log likelihood, as a function of the residual, is

\frac{\nu + 1}{2}\ln(1 + \lambda e_i^2/\nu)

plus a constant. For

\nu = 1

the t distribution is the Cauchy distribution and this is the Cauchy loss above, up to a constant factor, with

c^2 s^2 = 1/\lambda

. The scale mixture also explains the weights. Treat

\eta_i

as a latent variable. Its posterior given the residual is again a gamma distribution, with mean

\mathbb E[\eta_i \mid e_i] = \frac{\nu + 1}{\nu + \lambda e_i^2}.

An EM algorithm for the t regression computes these expectations in the E step and solves a weighted least squares problem with weights

\mathbb E[\eta_i \mid e_i]

in the M step. Bishop notes that the t fit can be found by EM (Section 2.3.7). The expectations are exactly the Cauchy weights for

\nu = 1

, so IRLS with Cauchy weights is EM for regression with Cauchy noise.

Comparison on contaminated data

The EM section modeled the same contamination differently: each target comes either from the line with Gaussian noise, with probability

\pi

, or from a uniform outlier distribution

\mathcal U

p(y^{(i)} \mid x^{(i)}) = \pi\, \mathcal N\big(y^{(i)} \mid \boldsymbol\theta^\top \boldsymbol\phi(x^{(i)}), \sigma^2\big) + (1 - \pi)\, \mathcal U(y^{(i)}).

The E step computes the responsibility

\gamma(z_i)

of the inlier component for each point, and the M step is a weighted least squares problem with weights

\gamma(z_i)

, followed by updates of

\sigma^2

and

\pi

, which has the same form as IRLS. The weight

\gamma(z_i)

falls toward zero much faster than a Cauchy weight: once a residual is a few

\sigma

from the line, the Gaussian density is negligible next to the uniform one. The mixture also returns two things a robust loss does not, an estimate of the inlier fraction

\pi

and a probability for each point that it is an inlier. The code fits all four methods to the 60-point data set of the EM section, in which 18 targets were replaced by gross errors.

def em_inlier_outlier(Phi, y, theta0, n_iter=200, pi0=0.5):
    # EM for the Gaussian plus uniform mixture; returns theta, sigma, pi and the responsibilities
    u = 1.0 / (y.max() - y.min())
    theta, sigma2, pi = theta0, np.var(y - Phi @ theta0), pi0
    for _ in range(n_iter):
        inlier = pi * norm.pdf(y, Phi @ theta, np.sqrt(sigma2))
        gamma = inlier / (inlier + (1 - pi) * u)
        theta = weighted_least_squares(Phi, y, gamma)
        sigma2 = max(np.sum(gamma * (y - Phi @ theta) ** 2) / gamma.sum(), 1e-12)
        pi = gamma.mean()
    return theta, np.sqrt(sigma2), pi, gamma

x_c, y_c, out_c = make_line_data(n_out=18)
Phi_c = design(x_c)
theta_ls_c = least_squares(Phi_c, y_c)
fits = {
    "least squares": theta_ls_c,
    "Huber": irls(Phi_c, y_c, huber_weight, theta_ls_c)[0],
    "Cauchy": irls(Phi_c, y_c, cauchy_weight, theta_ls_c)[0],
}
theta_em_c, sigma_em_c, pi_em_c, gamma_c = em_inlier_outlier(Phi_c, y_c, theta_ls_c)
fits["EM mixture"] = theta_em_c
fits["least squares, true inliers only"] = least_squares(Phi_c[~out_c], y_c[~out_c])
for name, th in fits.items():
    print(f"{name:34s} theta = [{th[0]:.2f}, {th[1]:.3f}]")
print(f"EM: sigma = {sigma_em_c:.2f}, pi = {pi_em_c:.2f}")

least squares                      theta = [2.10, 0.343]
Huber                              theta = [0.95, 0.506]
Cauchy                             theta = [0.87, 0.518]
EM mixture                         theta = [0.91, 0.511]
least squares, true inliers only   theta = [0.92, 0.508]
EM: sigma = 0.24, pi = 0.73

Scatter of 42 inliers on a line and 18 outliers spread over the plot, with four fitted lines. The least squares line is flatter than the others; the Huber, Cauchy and EM lines lie close together along the inliers.

Four fits to 60 points of which 18 are gross errors. Least squares is pulled to an intercept of 2.10 and a slope of 0.343. The three robust fits land within 0.01 of the slope of the least squares fit to the 42 true inliers, 0.508, which is the fit you would get if the outliers were known.

How much contamination each method tolerates

One data set says little about where each method fails. The code repeats the comparison on 50 random data sets for each outlier fraction from 0 to 0.7 and records the median absolute error of the slope. Each robust method starts from the least squares fit, which is what you would do without other information.

fractions = np.arange(0.0, 0.71, 0.05)
n_rep = 50
methods = ["least squares", "Huber", "Cauchy", "EM mixture"]
errors = {k: [] for k in methods}
for f in fractions:
    slopes = []
    for rep in range(n_rep):
        x, y, _ = make_line_data(n_out=int(round(f * 60)), seed=100 + rep)
        Phi = design(x)
        th_ls = least_squares(Phi, y)
        slopes.append([
            th_ls[1],
            irls(Phi, y, huber_weight, th_ls)[0][1],
            irls(Phi, y, cauchy_weight, th_ls)[0][1],
            em_inlier_outlier(Phi, y, th_ls)[0][1],
        ])
    for k, v in zip(methods, np.median(np.abs(np.array(slopes) - 0.5), axis=0)):
        errors[k].append(v)

print("outlier fraction", " ".join(f"{f:5.2f}" for f in fractions))
for k in methods:
    print(f"{k:16s}", " ".join(f"{v:5.3f}" for v in errors[k]))

outlier fraction  0.00  0.05  0.10  0.15  0.20  0.25  0.30  0.35  0.40  0.45  0.50  0.55  0.60  0.65  0.70
least squares    0.010 0.033 0.054 0.082 0.090 0.114 0.130 0.173 0.200 0.225 0.269 0.236 0.292 0.355 0.322
Huber            0.011 0.011 0.014 0.013 0.019 0.020 0.021 0.018 0.032 0.037 0.064 0.103 0.181 0.291 0.253
Cauchy           0.011 0.011 0.013 0.011 0.014 0.018 0.013 0.019 0.021 0.022 0.040 0.087 0.164 0.263 0.242
EM mixture       0.010 0.009 0.012 0.012 0.015 0.013 0.014 0.015 0.013 0.014 0.015 0.017 0.021 0.016 0.031

Median slope error against outlier fraction on a log scale for four methods. Least squares error grows steadily from the first outliers. Huber and Cauchy stay near the noise level up to about 0.45 and then rise sharply. The EM mixture stays near the noise level up to 0.7.

Median slope error over 50 data sets per outlier fraction, each robust method started from the least squares fit. Least squares degrades from the first outlier. Huber and Cauchy stay at the noise level, a median slope error between 0.01 and 0.04, up to an outlier fraction of 0.45, and then fail: at 0.6 their errors are 0.18 and 0.16, close to least squares. The scale estimate causes the failure. The median absolute deviation ignores at most half of the residuals. Beyond that, the outliers set

s

, every residual looks small relative to it, all weights approach 1, and IRLS returns to least squares. The EM mixture estimates

\sigma

from the points it believes are inliers and lets the uniform component absorb the rest, so its median error stays at 0.031 or below up to 0.7. These outliers are scattered at random, so they pull in no particular direction. Outliers that agree with each other are harder, as the next example shows.

Grid load in Winter Storm Uri

ERCOT operates the electricity grid that serves most of Texas. In mid-February 2021, Winter Storm Uri brought temperatures far below anything in the preceding winter, demand rose, and generation failed. From the early hours of 15 February, ERCOT shed load: it cut power to customers to keep the grid from collapsing. The load it reports for those hours is what it served, well below what customers would have drawn. A model of normal demand fitted to a window that includes those hours faces exactly the contamination of this section. The data are hourly ERCOT system load and a statewide temperature, from 1 December 2020 to 20 February 2021 (Central time). The load is the hourly native load from the ERCOT load archive. The temperature is the mean of the hourly readings at eight airports, one in each of ERCOT’s weather zones, from NOAA’s Integrated Surface Database. The table is published as the ercot-uri-load-temperature dataset, with the scripts that build it.

HF_REPO = "aegean-ai/ercot-uri-load-temperature"
grid = pd.read_parquet(f"https://huggingface.co/datasets/{HF_REPO}/resolve/main/data/ercot_hourly.parquet")
grid = grid[grid["window"].isin(["train", "control", "event"])].reset_index(drop=True)
local = grid["timestamp_local"]
ordinary = (grid["window"] == "train").to_numpy()  # 1 December 2020 to 7 February 2021
storm_week = (grid["window"] == "event").to_numpy()  # 14 to 19 February 2021
shed_days = local.dt.strftime("%m-%d").isin(["02-15", "02-16"]).to_numpy()
y_grid = grid["load_mw"].to_numpy() / 1000  # GW
print(f"{len(grid)} hours from {local.iloc[0]:%Y-%m-%d} to {local.iloc[-1]:%Y-%m-%d %H:%M}")
for name, mask in [("1 Dec to 7 Feb", ordinary), ("8 to 13 Feb", ~ordinary & ~storm_week), ("14 to 19 Feb", storm_week)]:
    print(f"{name:15s} temperature {grid['temp_f'][mask].min():5.1f} to {grid['temp_f'][mask].max():5.1f} F, "
          f"load {y_grid[mask].min():5.1f} to {y_grid[mask].max():5.1f} GW")

hours from 2020-12-01 to 2021-02-19 23:00
Dec to 7 Feb  temperature  28.4 to  77.5 F, load  30.9 to  58.6 GW
to 13 Feb     temperature  23.8 to  63.4 F, load  34.3 to  64.1 GW
to 19 Feb    temperature   6.2 to  42.9 F, load  39.8 to  69.7 GW

The model is the one used in the AI course’s grid assignment. Each hour is a row of a design matrix with 82 columns: the scaled temperature

T = (\text{temperature} - 60)/20

entered as

T

T^2

and

T^3

; 24 hour-of-day indicators and 7 day-of-week indicators; and the products of each hour indicator with

T

and with

T^2

, which give every hour of the day its own temperature response. The model is linear in

\boldsymbol\theta

, so every method of this section applies unchanged. The indicator columns and their products are not linearly independent (each set of indicators sums to one in every row), and the solver returns the minimum-norm solution, which does not change the fitted values.

def grid_design(df):
    # 82 columns: T, T^2, T^3, 24 hour and 7 day indicators, hour x T and hour x T^2
    T = (df["temp_f"].to_numpy() - 60.0) / 20.0
    H = np.eye(24)[df["hour"].to_numpy()]
    D = np.eye(7)[df["dow"].to_numpy()]
    return np.column_stack([T, T**2, T**3, H, D, H * T[:, None], H * (T**2)[:, None]])


Phi_grid = grid_design(grid)
print("design matrix:", Phi_grid.shape, " rank:", np.linalg.matrix_rank(Phi_grid))

design matrix: (1944, 82)  rank: 79

Two least squares fits serve as references. The first uses only the ordinary winter, 1 December to 7 February. The second uses the whole window, storm included. The code reports each fit’s root mean square residual on the ordinary winter, its mean residual on 15 and 16 February, and its mean residual on 14 February, the cold day before shedding began, when the load was still the demand.

def report(name, theta, gamma=None):
    e = y_grid - Phi_grid @ theta  # observed minus predicted, GW
    day14 = local.dt.strftime("%m-%d").eq("02-14").to_numpy()
    line = (f"{name:30s} RMS residual 1 Dec to 7 Feb {np.sqrt(np.mean(e[ordinary] ** 2)):4.2f} GW | "
            f"mean residual 15 to 16 Feb {e[shed_days].mean():+6.1f} GW | 14 Feb {e[day14].mean():+5.1f} GW")
    print(line)
    return e


theta_ordinary = least_squares(Phi_grid[ordinary], y_grid[ordinary])
theta_all = least_squares(Phi_grid, y_grid)
_ = report("least squares, ordinary winter", theta_ordinary)
_ = report("least squares, with the storm", theta_all)

least squares, ordinary winter RMS residual 1 Dec to 7 Feb 1.50 GW | mean residual 15 to 16 Feb  -39.2 GW | 14 Feb  -9.5 GW
least squares, with the storm  RMS residual 1 Dec to 7 Feb 1.74 GW | mean residual 15 to 16 Feb   -4.9 GW | 14 Feb  +9.7 GW

Fitted to the ordinary winter, the model predicts far more load on 15 and 16 February than ERCOT served: those residuals are the shed load. Fitted with the storm included, least squares bends to meet the shed hours, which shrinks their residuals and moves the error elsewhere. On 14 February, when the load was real demand, the bent model now under-predicts. The robust fits use the whole window. Huber and Cauchy IRLS start from the least squares fit of the whole window. EM is started twice: from the same fit, and from the fit to the ordinary winter.

theta_huber, _ = irls(Phi_grid, y_grid, huber_weight, theta_all)
theta_cauchy, _ = irls(Phi_grid, y_grid, cauchy_weight, theta_all)
theta_em_all, s_all, pi_all, gamma_all = em_inlier_outlier(Phi_grid, y_grid, theta_all, pi0=0.9)
theta_em_ord, s_ord, pi_ord, gamma_ord = em_inlier_outlier(Phi_grid, y_grid, theta_ordinary, pi0=0.9)
_ = report("Huber IRLS", theta_huber)
_ = report("Cauchy IRLS", theta_cauchy)
_ = report("EM, least squares start", theta_em_all)
_ = report("EM, ordinary-winter start", theta_em_ord)
print(f"EM, least squares start:   sigma = {s_all:.2f} GW, pi = {pi_all:.3f}")
print(f"EM, ordinary-winter start: sigma = {s_ord:.2f} GW, pi = {pi_ord:.3f}")

Huber IRLS                     RMS residual 1 Dec to 7 Feb 1.64 GW | mean residual 15 to 16 Feb   -3.8 GW | 14 Feb +11.9 GW
Cauchy IRLS                    RMS residual 1 Dec to 7 Feb 1.65 GW | mean residual 15 to 16 Feb   -2.5 GW | 14 Feb +13.3 GW
EM, least squares start        RMS residual 1 Dec to 7 Feb 1.63 GW | mean residual 15 to 16 Feb   -4.8 GW | 14 Feb +11.8 GW
EM, ordinary-winter start      RMS residual 1 Dec to 7 Feb 1.58 GW | mean residual 15 to 16 Feb  -23.7 GW | 14 Feb  -0.3 GW
EM, least squares start:   sigma = 1.52 GW, pi = 0.903
EM, ordinary-winter start: sigma = 1.48 GW, pi = 0.910

Hourly load from 11 to 19 February 2021 with four model curves. The served load drops sharply on 15 February. The least squares and Huber fits that include the storm dip toward the served load; the ordinary-winter fit and the EM fit stay well above it for most of 15 and 16 February, with brief dips in the EM curve.

Served load and fitted load from 11 to 19 February 2021. The shaded band marks 15 and 16 February.

day = local.dt.strftime("%m-%d")
late = (local >= "2021-02-08").to_numpy()
summary = pd.DataFrame({"day": day[late], "gamma": gamma_ord[late]}).groupby("day")["gamma"].agg(
    mean_responsibility="mean", hours_flagged=lambda g: int((g < 0.5).sum())
)
print(summary.round(3).to_string())
print(f"hours flagged in 1 Dec to 7 Feb: {(gamma_ord[ordinary] < 0.5).sum()} of {ordinary.sum()}")

       mean_responsibility  hours_flagged
day                                      
02-08                0.987              0
02-09                0.980              0
02-10                0.984              0
02-11                0.958              0
02-12                0.844              3
02-13                0.836              4
02-14                0.852              3
02-15                0.163             20
02-16                0.000             24
02-17                0.000             24
02-18                0.120             21
02-19                0.171             19
hours flagged in 1 Dec to 7 Feb: 20 of 1656

Inlier responsibility for each hour from 8 to 19 February 2021 for two EM runs. The run started from the ordinary winter drops to zero through 15, 16 and 17 February. The run started from least squares stays high through most of those hours.

Inlier responsibility $\gamma(z_i)$ per hour, 8 to 19 February, for the two EM starts. The shaded band marks 15 and 16 February.

h_grid = np.einsum("ij,ji->i", Phi_grid, np.linalg.pinv(Phi_grid))  # leverage h_ii
print(f"mean leverage, 1 Dec to 7 Feb: {h_grid[ordinary].mean():.3f}")
print(f"mean leverage, 15 to 16 Feb:   {h_grid[shed_days].mean():.3f}")
print(f"coldest hour before 8 Feb: {grid['temp_f'][ordinary].min():.1f} F; hours below it from 14 to 19 Feb: "
      f"{(grid['temp_f'][storm_week] < grid['temp_f'][ordinary].min()).sum()}")

mean leverage, 1 Dec to 7 Feb: 0.031
mean leverage, 15 to 16 Feb:   0.289
coldest hour before 8 Feb: 28.4 F; hours below it from 14 to 19 Feb: 126

Leverage and the starting point

Least squares fitted with the storm moves the mean residual of 15 and 16 February from

-39.2

GW to

-4.9

GW: the fit has bent to meet the shed hours. On 14 February, the cold day before shedding, the bent model now predicts 9.7 GW less than the real demand. Huber and Cauchy IRLS do not undo the bend. Their mean residuals on 15 and 16 February are

-3.8

and

-2.5

GW, because of leverage. The coldest hour of the ordinary winter was 28.4 F, and 126 hours of the storm week were colder. Those hours sit at the far end of the temperature axis, where the cubic and the per-hour quadratic terms are constrained by little else, and their mean leverage is 0.289 against 0.031 for the ordinary winter. The fit can pass close to them, and once it does their residuals are small and no weight function reduces their weight. It is the same effect as the single high-leverage point of the first experiment, produced here by a group of hours that agree with each other. The Huber loss is convex, so a different start would not change its answer. EM started from the least squares fit stops in the same place, a local maximum of the likelihood in which the storm hours are inliers. Started from the fit to the ordinary winter, EM keeps the shape of normal demand: its residual on 14 February is

-0.3

GW, and its residuals on 15 and 16 February stay large and negative. Its inlier responsibilities flag 20 of the 24 hours of 15 February, when shedding began in the early morning, and every hour of 16 and 17 February. Most hours of 18 and 19 February stay flagged as well (21 and 19 of 24). They also flag 7 hours of the cold days of 12 and 13 February. On all of these days the temperature fell below the range of the ordinary winter, where the model is extrapolating. The two EM runs use the same data and the same model and differ only in where they start. A robust method downweights points that disagree with its current fit. It cannot recover once the fit already agrees with the bad points, so it needs a start that is not already contaminated. Here that start came from knowing which window was ordinary. For models with few parameters, the next section, RANSAC, finds such a start with no prior knowledge, by fitting many small random subsets of the data. It also explains why that does not work for an 82-parameter model like this one. The example has three limits. The temperature is an unweighted mean over eight airports, so it does not weight the weather by where the load is. The load during shedding is what ERCOT served, so the size of the negative residuals estimates the unserved demand only as far as the model can be trusted at temperatures it was never fitted to. The results also come from a single storm. Key references: (Bishop, n.d.; Dempster et al., n.d.)

References

Dempster, A., Laird, N., Rdin, D. (n.d.). Maximum Likelihood from Incomplete Data via the EM Algorithm.
Bishop (n.d.). Pattern Recognition and Machine Learning.

Found an error on this page? Report an issue.

Foundations

Learning & Regression

Optimization

Maximum Likelihood

Classification

Dimensionality Reduction

Mixture Models and EM

Robust Estimation

Least Squares with Outliers

The running example

One bad point and leverage

Robust losses and iteratively reweighted least squares

The Student’s t view

Comparison on contaminated data

How much contamination each method tolerates

Grid load in Winter Storm Uri

Leverage and the starting point

References

​The running example

​One bad point and leverage

​Robust losses and iteratively reweighted least squares

​The Student’s t view

​Comparison on contaminated data

​How much contamination each method tolerates

​Grid load in Winter Storm Uri

​Leverage and the starting point

​References

The running example

One bad point and leverage

Robust losses and iteratively reweighted least squares

The Student’s t view

Comparison on contaminated data

How much contamination each method tolerates

Grid load in Winter Storm Uri

Leverage and the starting point

References