> ## Documentation Index > Fetch the complete documentation index at: https://aegean.ai/llms.txt > Use this file to discover all available pages before exploring further. # Least Squares with Outliers > How bad points bend a least squares fit, and how to resist them. 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](https://aegean.ai/aiml-common/lectures/mixture-of-gaussians/em-algorithm.md) 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](https://aegean.ai/aiml-common/lectures/regression/linear-regression.md) 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$. ```python theme={null} 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](https://aegean.ai/aiml-common/lectures/optimization/maximum-likelihood/conditional_maximum_likelihood.md) section shows. ```python theme={null} 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)) ``` ```output theme={null} 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$. ```python theme={null} 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") ``` ```output theme={null} 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.](https://mintcdn.com/aegeanaiinc/uADPagtO6esItILu/aiml-common/lectures/robust-estimation/least-squares-outliers/images/cell_5_output_1.svg?fit=max\&auto=format\&n=uADPagtO6esItILu\&q=85\&s=16696353197ce0e6d79b6e5f11efaada) *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](https://aegean.ai/aiml-common/lectures/rse/scan-registration.md) 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. ```python theme={null} 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$. ```python theme={null} 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}") ``` ```output theme={null} 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.](https://mintcdn.com/aegeanaiinc/uADPagtO6esItILu/aiml-common/lectures/robust-estimation/least-squares-outliers/images/cell_8_output_1.svg?fit=max\&auto=format\&n=uADPagtO6esItILu\&q=85\&s=e023228f0b57506fe58b436e34d08d73) *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. ```python theme={null} 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 ``` ```python theme={null} 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}") ``` ```output theme={null} 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.](https://mintcdn.com/aegeanaiinc/uADPagtO6esItILu/aiml-common/lectures/robust-estimation/least-squares-outliers/images/cell_11_output_1.svg?fit=max\&auto=format\&n=uADPagtO6esItILu\&q=85\&s=b5fc1f49875f2066f4e55f8915d93809) *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. ```python theme={null} 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])) ``` ```output theme={null} 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.](https://mintcdn.com/aegeanaiinc/uADPagtO6esItILu/aiml-common/lectures/robust-estimation/least-squares-outliers/images/cell_13_output_1.svg?fit=max\&auto=format\&n=uADPagtO6esItILu\&q=85\&s=f191a8885fb35684d51e626411b4f837) *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](https://huggingface.co/datasets/aegean-ai/ercot-uri-load-temperature) dataset, with the scripts that build it. ```python theme={null} 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") ``` ```output theme={null} 1944 hours from 2020-12-01 to 2021-02-19 23:00 1 Dec to 7 Feb temperature 28.4 to 77.5 F, load 30.9 to 58.6 GW 8 to 13 Feb temperature 23.8 to 63.4 F, load 34.3 to 64.1 GW 14 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. ```python theme={null} 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)) ``` ```output theme={null} 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. ```python theme={null} 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) ``` ```output theme={null} 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. ```python theme={null} 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}") ``` ```output theme={null} 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.](https://mintcdn.com/aegeanaiinc/uADPagtO6esItILu/aiml-common/lectures/robust-estimation/least-squares-outliers/images/cell_18_output_1.svg?fit=max\&auto=format\&n=uADPagtO6esItILu\&q=85\&s=33979e7ba90b5015379c3c87bb7d7a7a) *Served load and fitted load from 11 to 19 February 2021. The shaded band marks 15 and 16 February.* ```python theme={null} 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()}") ``` ```output theme={null} 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.](https://mintcdn.com/aegeanaiinc/uADPagtO6esItILu/aiml-common/lectures/robust-estimation/least-squares-outliers/images/cell_20_output_1.svg?fit=max\&auto=format\&n=uADPagtO6esItILu\&q=85\&s=e469cdcbd9c79376705fa86d53dd7512) *Inlier responsibility $\gamma(z_i)$ per hour, 8 to 19 February, for the two EM starts. The shaded band marks 15 and 16 February.* ```python theme={null} 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()}") ``` ```output theme={null} 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](https://aegean.ai/aiml-common/lectures/robust-estimation/ransac.md), 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](https://github.com/aegean-ai/eaia-feedback/issues/new?template=page-issue.yml\&page=/aiml-common/lectures/robust-estimation/least-squares-outliers). This documentation is built and hosted on [Mintlify](https://mintlify.com), a developer documentation platform.