Loss Functions - What You Are Actually Minimizing

Context

A loss function is the formula behind the one number a fit drives down, and the same model fitted to the same listings under two losses comes out with two different sets of parameters.

You meet one every time a price model is fitted: deciding whether a model should chase the listings priced far above the rest, reading an estimator's documentation to learn which center its line aims at, and stopping a teammate from pasting a squared price into a report as if it were a miss in currency.

The failures are quiet. Two fits of one column return two slopes and the difference passes for noise, a mistyped price drags one model and leaves another exactly where it was, a loss of 38 billion is read as currency, and a loss computed on other rows announces a change no fit made. Lesson 3.1.2 scored a fit on listings it never saw, and this lesson opens the fit itself to name what it was minimizing all along.

Goal

• Scan a model's parameter over a grid of candidates with NumPy, find the smallest loss with np.argmin, and check the answer against np.linalg.lstsq.
• Explain precisely, with no derivative, why squared error chooses the mean price and absolute error the median.
• Fit one feature under both losses with LinearRegression and QuantileRegressor, and trace how far each line moves when the highest prices change.
• Draw a loss curve on a fig, ax figure, save it, and prove what it holds by printing it.
• Distinguish the loss a fit minimizes from the number you report, and report a fit by its average miss in currency on held-out rows.

Prerequisites

None of this is required, and the lesson stands on its own. Both are here if you want the idea drawn before you write it in NumPy.

• YouTube - "The Main Ideas of Fitting a Line to Data (The Main Ideas of Least Squares and Linear Regression.)" by StatQuest with Josh Starmer. A line fitted to a handful of points by squaring each residual, adding the squares and keeping the line whose sum is smallest.
• Google Machine Learning Crash Course - "Linear regression: Loss". Absolute and squared loss defined side by side, and a "Choosing a loss" section on how differently each one treats a point far from the rest.

Table of Contents

• What a Loss Function Is
• Squared Error, and Why the Mean Wins It
• Absolute Error, and Why the Median Wins It
• The Loss Decides the Fit
• The Loss You Minimize and the Number You Report

Theory

Every block in this section adds to one script, losses.py in ~/dev/bootcamp, run with python3 losses.py beside listings.csv.

What a Loss Function Is

What. The loss Lesson 3.1.1 named is a function of the parameters. With the rows held fixed, it turns each candidate's predictions into one number, and the fit looks for the candidate at its lowest point. Why. Every choice a fit makes follows from that number, so whatever the loss rewards is what the fitted parameters deliver, and two losses on the same rows can deliver two different models. When. Every time a model is fitted, whether you chose the loss or the estimator chose it for you.

The model with the fewest moving parts is a price per square meter and nothing else, price = w * area, 2.3.4's line through the origin, now on listings.csv. Its loss here is the mean squared error, each listing's residual squared and averaged over the listings, and since only w moves while the search runs, the first block writes the loss as a function of w.

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd

listings = pd.read_csv("listings.csv")
# the 125 listings with no area dropped, as in 3.1.1
has_area = listings.dropna(subset=["area_sqm"])
area = has_area["area_sqm"].to_numpy()
price = has_area["price"].to_numpy()
print(len(price))            # 2875


def mse(w):
    predicted = w * area     # the model: price = w * area
    return np.mean((price - predicted) ** 2)


print(f"{mse(3_000):,.0f}")  # 24,320,084,049
print(f"{mse(4_000):,.0f}")  # 18,456,928,403

Finding the best w takes a scan (trying every candidate on a grid, an evenly spaced range of values, and keeping the one with the smallest loss), and 2.3.3's np.argmin returns the position of the smallest.

grid = np.arange(0, 6_001, 1.0)  # 0, 1, 2, ..., 6000
losses = np.array([mse(w) for w in grid])
best = grid[np.argmin(losses)]   # the candidate with the smallest loss
print(best)                      # 3750.0
print(f"{losses.min():,.0f}")    # 17,720,312,891

# 2.3.4's projection through the origin
w_ls, ssr, rank, sv = np.linalg.lstsq(area[:, np.newaxis], price)
print(np.round(w_ls, 2))         # [3749.58]

The grid holds whole numbers, so 3,750 is the best whole number, and lstsq reaches the exact bottom, 3,749.58, in one step where the scan tried 6,001 candidates. 2.3.4 minimized the sum of squared residuals instead of their mean, and dividing every candidate's sum by the same 2,875 leaves the bottom where it was. A scan works for one parameter or two and is hopeless for thousands, which is why Module 4 searches without a grid.

The curve is drawn as 2.7.2 drew its charts, on a fig, ax figure saved with savefig, with the ax.axvline of 2.7.2's Assignment at the best candidate, and checked by printing what it holds:

fig, ax = plt.subplots(figsize=(8, 5))
ax.plot(grid, losses)
ax.axvline(best, linestyle="--", color="black", label=f"lowest loss at w = {best:.0f}")
ax.set_xlabel("w, price per square meter")
ax.set_ylabel("mean squared error")
ax.legend()
fig.savefig("loss_curve.png", dpi=150, bbox_inches="tight")

curve = ax.lines[0]                    # the loss curve, drawn first
drawn_w, drawn_loss = curve.get_xdata(), curve.get_ydata()
print(len(ax.lines))                   # 2
print(len(drawn_w))                    # 6001
print(drawn_w[np.argmin(drawn_loss)])  # 3750.0
plt.close(fig)

The Axes holds the curve and the marker, and the curve drew all 6,001 candidates, falling to one lowest point at 3,750 and rising on both sides of it.

DON'T:

• DON'T trust a scan whose answer sits on the edge of its grid. A scan over np.arange(0, 3_001, 1.0) returns 3000.0, its own last candidate, with no warning, because the loss is still falling where the grid stops. Widen the grid until the answer lands inside it.
• DON'T score a fit with raw residuals. np.mean(price - w * area) lets misses above and below cancel and falls without end as w grows, so the same scan returns its last candidate, 6000.0. A loss charges every miss, as squaring and absolute value both do.

Squared Error, and Why the Mean Wins It

What. Squared error charges each listing its residual squared, and the mean squared error averages the charges. For the simplest model, one price predicted for every listing, it chooses the mean price. Why. Module 2 fitted least squares without saying why squares. A miss twice as large costs four times as much, so the loss bends a fit toward the listings it misses most, and its lowest point is the one least squares reaches in one step, 2.3.4's projection, which is why LinearRegression minimizes it. When. Whenever a large miss costs far more than several small ones, and every time you fit a LinearRegression, which minimizes squared error whether you chose it or not.

A model that predicts one price for every listing needs no feature, so it uses all 3,000 listings, and it is the linear model whose design matrix is a column of ones and nothing else, so 2.3.4's projection reaches its answer without a single candidate:

all_prices = listings["price"].to_numpy()  # all 3,000, none missing
candidates = np.arange(0, 1_000_001, 1_000)
sq_losses = np.array([np.mean((all_prices - c) ** 2) for c in candidates])
print(candidates[np.argmin(sq_losses)])    # 381000
print(all_prices.mean())                   # 380763.0

ones = np.ones((len(all_prices), 1))       # one column, all ones
w_c, ssr, rank, sv = np.linalg.lstsq(ones, all_prices)
print(w_c)                                 # [380763.]

The scan's best is the candidate nearest the mean, and the projection lands on the mean itself. Two sets of four misses that add up to the same 200,000 show what squaring does:

spread_out = np.array([50_000, 50_000, 50_000, 50_000])
one_big = np.array([0, 0, 0, 200_000])
print(f"{np.mean(spread_out ** 2):,.0f}")  # 2,500,000,000
print(f"{np.mean(one_big ** 2):,.0f}")     # 10,000,000,000

# the mean's squared misses, largest first
misses = np.sort((all_prices - all_prices.mean()) ** 2)[::-1]
share = misses[:30].sum() / misses.sum()   # the 30 largest, 1 in 100
print(round(share, 3))                     # 0.177

One miss of 200,000 costs four times as much as four misses of 50,000, so a fit under squared error would rather miss many listings a little than one a lot, and on the listings, the 30 prices the mean misses most carry 17.7% of its loss.

Why the mean wins needs no derivative. Write each price's distance from a candidate c as its distance from the mean plus the mean's distance from c: squaring and averaging leaves the two squares and a middle term that averages to zero, because distances from the mean average to zero. The loss splits into np.var(all_prices), 2.2.5's mean squared distance from the mean, plus (c - all_prices.mean()) ** 2:

spread = np.var(all_prices)  # 2.2.5's mean squared distance from the mean
print(f"{spread:,.0f}")      # 38,129,487,498
for c in [300_000, 381_000, 450_000]:
    loss = np.mean((all_prices - c) ** 2)
    print(c, np.isclose(loss, spread + (c - all_prices.mean()) ** 2))

300000 True
381000 True
450000 True

The first part is the spread of the prices around their own mean, the same for every candidate, and no single price removes it. The second is zero at the mean and grows with the squared distance from it, so the mean wins and the curve rises alike on both sides, which is why the scan picked 381,000, 237 away, over 380,000, 763 away. The unit is squared currency, and Lesson 3.1.4 takes the root that turns a squared-error number back into currency.

DON'T:

• DON'T read a squared-error number in currency. 38,129,487,498 is a mean of squared misses, and no listing is off by anything near 38 billion.
• DON'T expect squared error to weigh every listing alike. The 30 listings the mean misses most carry 17.7% of its loss, and a fit under squared error leans toward them.

Absolute Error, and Why the Median Wins It

What. Absolute error charges each listing the size of its residual with the sign dropped, and the mean absolute error averages the charges, Lesson 3.1.2's average miss in currency used as a loss. The one price it chooses for every listing is the median. Why. Moving a prediction by a step changes the absolute charge of every listing the step does not cross by exactly that step, whatever the listing's distance, so only the side a listing sits on matters, and a price 1,500,000 above the rest weighs exactly as much as one 1,000 above. When. When a miss twice as large costs twice as much and no more, and when the target has a long tail, prices running far above the middle, that you do not want the fit to chase.

abs_losses = np.array([np.mean(np.abs(all_prices - c)) for c in candidates])
print(candidates[np.argmin(abs_losses)])  # 337000
print(np.median(all_prices))              # 337000.0

The reason can be counted. Every price in the file is a whole multiple of 1,000, so no listing sits strictly between 300,000 and 301,000, and a step up from 300,000 brings each listing above it 1,000 closer and takes each listing at or below it 1,000 farther, whatever its distance.

print(np.sum(all_prices % 1_000 != 0))  # 0
c = 300_000
at_or_below = np.sum(all_prices <= c)
above = np.sum(all_prices > c)
print(at_or_below, above)               # 1211 1789
here = np.mean(np.abs(all_prices - c))
up = np.mean(np.abs(all_prices - (c + 1_000)))
print(round(up - here, 2))              # -192.67
counted = 1_000 * (at_or_below - above) / len(all_prices)
print(round(counted, 2))                # -192.67
print(np.sum(all_prices < 337_000))     # 1493
print(np.sum(all_prices == 337_000))    # 9
print(np.sum(all_prices > 337_000))     # 1498

The loss fell by 192.67 because 578 more listings sit above 300,000 than at or below it, and the count predicts the step to the cent. It keeps falling while the side above outnumbers the side at or below, and stops where the sides balance. With 1,493 prices below 337,000, 9 on it and 1,498 above, a step up takes 1,502 listings farther and brings 1,498 closer, a step down takes 1,507 farther and brings 1,493 closer, and both directions cost: 337,000 is the bottom.

Squared error answers the same step differently. A listing whose residual is r has its squared charge changed by 1_000 ** 2 - 2 * 1_000 * r, which grows with its distance: from 300,000, the listing priced at 1,930,000 sheds 3,259,000,000 on that one step and a listing 1,000 above the candidate sheds 1,000,000, while under absolute error both shed exactly 1,000. That is the whole difference between the two centers, and the prices above the median reach much farther than the prices below it:

mid = np.median(all_prices)
p10, p90 = np.percentile(all_prices, [10, 90])
print(all_prices.max() - mid)  # 1593000.0
print(mid - all_prices.min())  # 277000.0
print(round(p90 - mid))        # 311200
print(round(mid - p10))        # 157000

The 90th percentile sits about twice as far above the median as the 10th sits below it. Squared error weighs every one of those distances and absolute error only counts the listings on each side, so that long upper tail holds the mean 43,763 above the median. A single listing shows the same effect:

typo = all_prices.copy()
typo[typo.argmax()] *= 10    # 1,930,000 keyed in as 19,300,000
mean_shift = typo.mean() - all_prices.mean()
median_shift = np.median(typo) - np.median(all_prices)
print(f"{mean_shift:,.0f}")  # 5,790
print(median_shift)          # 0.0

The mistyped listing sat above the median before the typo and sits above it after, which is all absolute error asks of it, while the mean takes a 3,000th of the extra 17,370,000. Both curves, side by side, each marked at its lowest candidate:

fig, axes = plt.subplots(1, 2, figsize=(14, 5))
for ax, curve, name in zip(axes, [sq_losses, abs_losses], ["squared", "absolute"]):
    lowest = candidates[np.argmin(curve)]
    ax.plot(candidates, curve)
    ax.axvline(lowest, linestyle="--", color="black", label=f"lowest at {lowest:,}")
    ax.set_xlabel("one price predicted for every listing")
    ax.set_ylabel(f"mean {name} error")
    ax.legend()
fig.savefig("constant_losses.png", dpi=150, bbox_inches="tight")

for ax in axes:
    print(len(ax.lines), ax.lines[1].get_xdata()[0])
plt.close(fig)

2 381000
2 337000

Each panel holds two lines, its curve and its marker, and the markers sit at the two centers.

DON'T:

• DON'T take the mean price as the price of a typical listing. 1,785 of the 3,000 listings sit below 380,763, and one mistyped price moved it by 5,790 while the median stayed where it was.
• DON'T count a listing that sits exactly on the candidate with the ones above it. 11 listings are priced at exactly 300,000, a step up takes them 1,000 farther, and a count made with < and >= predicts -200.0 instead of the real -192.67.
• DON'T read the median's stillness as absolute error ignoring the far listings. Each still counts once, on its side: drop the 30 highest prices and the median falls from 337,000 to 335,500, while the mean falls by 8,004.

The Loss Decides the Fit

What. The same model on the same feature, fitted under two losses, gives two lines. LinearRegression minimizes squared error, and scikit-learn's QuantileRegressor (a linear model that predicts a chosen quantile of the target, 2.2.5's cut point, rather than its mean) minimizes absolute error at quantile=0.5, which makes its line median regression. Why. A fitted line keeps its loss's center. A squared-error line with an intercept leaves residuals that average to zero, as the mean does, and an absolute-error line with an intercept leaves no more than half the listings on either side of it, as the median does. When. Every time you choose an estimator, because you choose its loss with it, and above all on a target with a long tail, such as price.

Two hyperparameters set QuantileRegressor up. quantile=0.5 asks for the median, and at that quantile the loss it minimizes is half the absolute error, which has the same lowest point. alpha=0.0 switches off a penalty on the size of the slope (an extra charge added to the loss that grows with the slope), whose strength is 1.0 unless you set it.

from sklearn.linear_model import LinearRegression, QuantileRegressor

X = has_area[["area_sqm"]]
y = has_area["price"]
squared_fit = LinearRegression().fit(X, y)
absolute_fit = QuantileRegressor(quantile=0.5, alpha=0.0).fit(X, y)
for fit in [squared_fit, absolute_fit]:
    residuals = y - fit.predict(X)
    sq_loss = np.mean(residuals ** 2)
    abs_loss = np.mean(np.abs(residuals))
    slope, intercept = round(fit.coef_[0], 1), round(fit.intercept_)
    print(slope, intercept, f"{sq_loss:,.0f}", f"{abs_loss:,.0f}")

3288.3 54717 17,225,527,337 92,547
3146.1 45596 17,803,356,831 90,751

On the same 2,875 listings each line wins its own loss and loses the other's: LinearRegression's has the smaller mean squared error, QuantileRegressor's the smaller mean absolute error. Neither is the better line, each is the best line under its loss.

residuals = y - absolute_fit.predict(X)
# above the median line, below it, and on it to within 1
print(np.sum(residuals > 1))           # 1436
print(np.sum(residuals < -1))          # 1437
print(np.sum(np.abs(residuals) <= 1))  # 2
print(round(residuals.mean()))         # 23207
higher = squared_fit.predict(X) > absolute_fit.predict(X)
print(higher.sum())                    # 2875

The median line passes through 2 listings and leaves 1,436 above it and 1,437 below, the median's balance of counts carried along the area. The listings above reach farther than the ones below, so its residuals average 23,207, while the squared-error line's average zero, 2.3.4's column of ones at work, and the squared-error line runs above the median line at the area of every one of the 2,875 listings.

Which line is steeper is a separate question, and the far listings are exactly what the two losses weigh differently. Three of them show how differently:

# the 2,872 listings left without the three highest prices
trimmed = has_area.sort_values("price").iloc[:-3]
X_t, y_t = trimmed[["area_sqm"]], trimmed["price"]
squared_t = LinearRegression().fit(X_t, y_t)
absolute_t = QuantileRegressor(quantile=0.5, alpha=0.0).fit(X_t, y_t)
print(round(squared_t.coef_[0], 1))   # 3249.2
print(round(absolute_t.coef_[0], 1))  # 3141.8

Three listings out of 2,875 moved the squared-error slope by 39.1 and the median line's by 4.3. Squared error charged each of the three its distance from the line squared, and to the median line they were 3 of the 1,436 listings above it.

QuantileRegressor cannot reach its line by a projection, because absolute error has no one-step formula like 2.3.4's, so it hands the problem to scipy.optimize.linprog, SciPy's linear-programming solver (a routine that finds the best values for a problem written as weighted sums under limits). That routine is one reason Lesson 3.1.1 installed SciPy.

DON'T:

• DON'T assume every estimator minimizes the same loss. On the same 2,875 listings LinearRegression and QuantileRegressor returned slopes of 3,288.3 and 3,146.1, each the best line under its own loss.
• DON'T leave QuantileRegressor's alpha at its default. QuantileRegressor(quantile=0.5) returns a slope of 2,966.0, its penalty of 1.0 pulling the slope down, and its mean absolute error of 90,942 sits above the 90,751 of alpha=0.0, so it is not the median line.
• DON'T compare losses computed on different rows. squared_fit, unchanged, scores a mean squared error of 17,225,527,337 on all 2,875 listings and 15,999,788,517 on the 2,872 in trimmed: the line did not improve, three hard listings left.

The Loss You Minimize and the Number You Report

What. The loss is the number a fit drives down on its training rows. The number you report is computed afterward, on held-out rows, to answer the question of the person who reads it, and a fit is routinely reported by a number its loss never was. Why. A loss is chosen for how it behaves in the search, a single bottom and a reliable way down to it. A reported number is chosen for what its reader will act on, and a pricing team asks how far off a typical price is, in currency, on listings the model has not seen. When. Every time a result leaves your script, in a report, a dashboard or a message to the pricing team.

Both models are fitted on the training rows of Lesson 3.1.2's split and reported by 3.1.2's mean absolute error on the test rows:

from sklearn.model_selection import train_test_split

X_train, X_test, y_train, y_test = train_test_split(
    X, y, test_size=0.2, random_state=42
)
model = LinearRegression().fit(X_train, y_train)
median_model = QuantileRegressor(quantile=0.5, alpha=0.0).fit(X_train, y_train)

train_loss = np.mean((y_train - model.predict(X_train)) ** 2)
print(f"{train_loss:,.0f}")  # 18,127,720,025
for fitted in [model, median_model]:
    test_mae = np.abs(y_test - fitted.predict(X_test)).mean()
    print(f"{test_mae:,.0f}")

87,952
85,031

The first number is what LinearRegression minimized, in squared currency, on the 2,300 listings it was fitted to. It never computed an average miss, and the pricing team never asked for squared currency. The other two answer the team on 575 listings neither model saw: a price is off by 87,952 on average under the squared-error fit and by 85,031 under the median line, the fit whose loss matches the question. One split is one draw, so that gap belongs to this split.

What keeps squared error the usual loss, even when the report is an average miss, is how it behaves in the search. The 2,300 training prices are an even count, and their two middle prices differ:

train_prices = np.sort(y_train.to_numpy())
middle = train_prices[1149:1151]   # the 1,150th and 1,151st
print(middle)                      # [336000 337000]
flat = [np.mean(np.abs(train_prices - c)) for c in [336_000, 336_500, 337_000]]
print(np.allclose(flat, flat[0]))  # True

Every constant from 336,000 to 337,000 has the same absolute error on these rows: the bottom is a flat stretch 1,000 wide, the loss names no single best price, and which point a search returns depends on the search. Squared error's bottom for one price is always the single point at the mean, reached in one step by a projection. Those two properties are what squared error is chosen for, and neither says anything about what the pricing team wants to read, which is why the fit and the report may use different numbers.

DON'T:

• DON'T report a training loss. 18,127,720,025 is squared currency measured on the listings the fit was chosen on, and it answers no question the pricing team asked.
• DON'T expect np.argmin to tell you a bottom is flat. Scanned on the grid of 1,000, the training prices' absolute error ties at 336,000 and 337,000, and np.argmin returns the first, 336,000, without a word.

Lab Exercises

Lab Exercise 1: A House by the Square Meter, and by One Price

Work in house_loss.py in ~/dev/bootcamp, on the houses in listings.csv that have an area.

Scan the price per square meter w of price = w * area from 0 to 6,000 in steps of 1 under mean squared error, and print the best w, its loss and the exact answer from np.linalg.lstsq through the origin, to two decimals. Draw the loss curve with a dashed vertical line at the best w, save it as house_loss.png, and print how many lines the Axes holds, how many points the curve drew and the w at its lowest drawn point.

Then scan one price for every house, 0 to 1,500,000 in steps of 1,000, under mean squared error. Print the best candidate, the houses' mean price and the constant lstsq finds on a column of ones, and confirm that the loss at the best candidate equals the houses' variance plus the squared distance from the candidate to their mean. Finally print the share of the constant's loss carried by the 1 in 100 houses the mean misses most, and the constant's loss itself, with a comment comparing it with the line's.

Lab Exercise 2: Eastgate Under Two Losses

Work in eastgate_loss.py in ~/dev/bootcamp, on the Eastgate listings that have an area.

Scan one price for every listing, 0 to 2,000,000 in steps of 1,000, under mean absolute error, and print the best candidate and the median. At 500,000, count the listings at or below the candidate and above it, and show that moving the candidate to 501,000 changes the loss by what the counts predict.

Fit LinearRegression and QuantileRegressor(quantile=0.5, alpha=0.0) on area_sqm, and print each fit's slope, intercept, mean squared error and mean absolute error. Drop the three highest prices, refit both, and print each new slope and how far it moved. Finally split the listings with test_size=0.2 and random_state=42, fit both models on the training rows, print the squared-error model's training loss and each model's mean absolute error on the test rows, and add a comment saying which numbers go to the pricing team.

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd

listings = pd.read_csv("listings.csv")
houses = listings[listings["property_type"] == "house"]
houses = houses.dropna(subset=["area_sqm"])
area = houses["area_sqm"].to_numpy()
price = houses["price"].to_numpy()
print(len(price))              # 770

grid = np.arange(0, 6_001, 1.0)
losses = np.array([np.mean((price - w * area) ** 2) for w in grid])
best = grid[np.argmin(losses)]
w_ls, ssr, rank, sv = np.linalg.lstsq(area[:, np.newaxis], price)
print(best)                    # 3515.0
print(f"{losses.min():,.0f}")  # 27,999,084,702
print(np.round(w_ls, 2))       # [3514.56]

fig, ax = plt.subplots(figsize=(8, 5))
ax.plot(grid, losses)
ax.axvline(best, linestyle="--", color="black", label=f"lowest loss at w = {best:.0f}")
ax.set_xlabel("w, price per square meter")
ax.set_ylabel("mean squared error")
ax.legend()
fig.savefig("house_loss.png", dpi=150, bbox_inches="tight")

drawn_w, drawn_loss = ax.lines[0].get_xdata(), ax.lines[0].get_ydata()
print(len(ax.lines), len(drawn_w))     # 2 6001
print(drawn_w[np.argmin(drawn_loss)])  # 3515.0
plt.close(fig)

candidates = np.arange(0, 1_500_001, 1_000)
sq_losses = np.array([np.mean((price - c) ** 2) for c in candidates])
best_c = candidates[np.argmin(sq_losses)]
ones = np.ones((len(price), 1))
w_c, ssr, rank, sv = np.linalg.lstsq(ones, price)
print(best_c)                     # 528000
print(round(price.mean(), 2))     # 527529.87
print(np.round(w_c, 2))           # [527529.87]
distance = (best_c - price.mean()) ** 2
# True: the houses' variance plus the squared distance
print(np.isclose(sq_losses.min(), np.var(price) + distance))

misses = np.sort((price - price.mean()) ** 2)[::-1]
top = len(price) // 100           # 7 houses, 1 in 100
share = misses[:top].sum() / misses.sum()
print(round(share, 3))            # 0.126
print(f"{sq_losses.min():,.0f}")  # 46,584,851,948
# the line's 27,999,084,702 on the same 770 houses is far lower

import numpy as np
import pandas as pd
from sklearn.linear_model import LinearRegression, QuantileRegressor
from sklearn.model_selection import train_test_split

listings = pd.read_csv("listings.csv")
eastgate = listings[listings["neighborhood"] == "Eastgate"]
eastgate = eastgate.dropna(subset=["area_sqm"])
prices = eastgate["price"].to_numpy()
print(len(prices))                        # 331

candidates = np.arange(0, 2_000_001, 1_000)
abs_losses = np.array([np.mean(np.abs(prices - c)) for c in candidates])
print(candidates[np.argmin(abs_losses)])  # 542000
print(np.median(prices))                  # 542000.0

c = 500_000
at_or_below = np.sum(prices <= c)
above = np.sum(prices > c)
here = np.mean(np.abs(prices - c))
up = np.mean(np.abs(prices - (c + 1_000)))
print(at_or_below, above)                 # 136 195
print(round(up - here, 2))                # -178.25
counted = 1_000 * (at_or_below - above) / len(prices)
print(round(counted, 2))                  # -178.25

X, y = eastgate[["area_sqm"]], eastgate["price"]
squared_fit = LinearRegression().fit(X, y)
absolute_fit = QuantileRegressor(quantile=0.5, alpha=0.0).fit(X, y)
for fit in [squared_fit, absolute_fit]:
    residuals = y - fit.predict(X)
    slope, intercept = round(fit.coef_[0], 1), round(fit.intercept_)
    sq_loss = np.mean(residuals ** 2)
    abs_loss = np.mean(np.abs(residuals))
    print(slope, intercept, f"{sq_loss:,.0f}", f"{abs_loss:,.0f}")

trimmed = eastgate.sort_values("price").iloc[:-3]
X_t, y_t = trimmed[["area_sqm"]], trimmed["price"]
squared_t = LinearRegression().fit(X_t, y_t)
absolute_t = QuantileRegressor(quantile=0.5, alpha=0.0).fit(X_t, y_t)
for before, after in [(squared_fit, squared_t), (absolute_fit, absolute_t)]:
    moved = after.coef_[0] - before.coef_[0]
    print(round(after.coef_[0], 1), round(moved, 1))

4728.0 78278 20,548,697,694 87,099
4803.0 53252 20,857,415,476 85,308
4511.1 -216.9
4688.7 -114.3

X_train, X_test, y_train, y_test = train_test_split(
    X, y, test_size=0.2, random_state=42
)
squared_model = LinearRegression().fit(X_train, y_train)
median_model = QuantileRegressor(quantile=0.5, alpha=0.0).fit(X_train, y_train)
train_loss = np.mean((y_train - squared_model.predict(X_train)) ** 2)
print(f"{train_loss:,.0f}")
for model in [squared_model, median_model]:
    test_mae = np.abs(y_test - model.predict(X_test)).mean()
    print(f"{test_mae:,.0f}")
# the two test errors go to the pricing team, in currency, on 67
# listings neither model saw, and the training loss stays here

21,968,958,631
84,189
80,300

Theory Assessment

Q1: Three Answers to One Scan

The scan of w returned 3,750.0, np.linalg.lstsq returned 3,749.58, and a teammate's scan over np.arange(0, 3_001, 1.0) returned 3,000.0. Explain each number, and say whether scanning the sum of the squared residuals, 2.3.4's quantity, instead of their mean would change any of them.

Hint: Ask which values a scan is able to return, and what the lowest candidate is on a curve still falling where the grid ends. For the sum, ask what dividing every candidate's loss by the same 2,875 does to their order.

Q2: The Two Parts of a Squared Loss

Predicting 381,000 for every listing scores a mean squared error of 38,129,543,667 on the 3,000 prices, and predicting their mean scores 38,129,487,498. Account for the difference of 56,169 with the split of the squared loss into two parts, say why no single price can score below 38,129,487,498, and say what that number measures and in what unit.

Hint: One part of the split changes with the candidate and the other does not. Start from 381,000 minus 380,763.

Q3: Counting Sides

Moving the constant from 300,000 to 301,000 lowered the mean absolute error by 192.67. Explain that number from the 1,211 listings at or below 300,000 and the 1,789 above it. Then use the same reasoning to explain why a mistyped price moved the mean by 5,790 and the median by nothing, and why dropping the three highest prices moved LinearRegression's slope by 39.1 and QuantileRegressor's by 4.3.

Hint: Under absolute error, ask how much one listing's charge changes when the prediction moves 1,000, and whether the listing's distance enters. Under squared error, ask what 1_000 ** 2 - 2 * 1_000 * r does as r grows.

Q4: The Typo That Moved One Line

A data-entry error multiplies the highest price in the table, 1,930,000, by 10. Read this script without running it.

import numpy as np
import pandas as pd
from sklearn.linear_model import LinearRegression, QuantileRegressor

listings = pd.read_csv("listings.csv")
clean = listings.dropna(subset=["area_sqm"]).sort_values("price")
typo = clean.copy()
top = typo.index[-1]  # the label of the highest price
typo.loc[top, "price"] = typo.loc[top, "price"] * 10

for frame in [clean, typo]:
    X, y = frame[["area_sqm"]], frame["price"]
    squared_line = LinearRegression().fit(X, y)
    median_line = QuantileRegressor(quantile=0.5, alpha=0.0).fit(X, y)
    median_loss = np.mean(np.abs(y - median_line.predict(X)))
    print(round(squared_line.coef_[0], 1), round(median_line.coef_[0], 1), round(median_loss))

Say what the two lines print. A teammate then argues that the median line must have changed because its loss rose by about 6,042, and that the squared-error fit is ruined because its mean squared error rose from 17,225,527,337 to 140,134,945,403. Answer both arguments.

Hint: For every line anywhere near the best, ask how much the typo raises the mistyped listing's absolute charge, and whether that amount depends on which line it is. Then ask what rows the two mean squared errors were computed on.

3288.3 3146.1 90751
3344.6 3146.1 96793

Assignment

Assignment: Huber Loss, Between the Mean and the Median

The pricing team wants a loss that charges an ordinary miss like squared error and the long tail's misses like absolute error, so a few far listings stop dragging the fit while small misses still cost more than in proportion. Huber loss does this with a threshold delta, in currency: a residual r no larger than delta in size is charged r ** 2 / 2, and a larger one delta * (abs(r) - delta / 2), a straight line that meets the curve at the threshold with no jump.

In huber_loss.py in ~/dev/bootcamp, write huber(residuals, delta) with NumPy, returning one charge per residual. On all 3,000 prices, scan the lesson's candidates for one price predicted for every listing under the mean Huber loss with delta at 1,000, 50,000, 150,000, 500,000 and 2,000,000, and print each best candidate and how many listings sit past the threshold at it, with the median and the mean beside them. Prove that at 2,000,000 the mean Huber loss of every candidate is exactly half its mean squared error, and explain from the two formulas why the scan at 1,000 lands where absolute error did. Then scan w of price = w * area on the 2,875 listings with an area, from 0 to 6,000 in steps of 1, under squared error, absolute error and Huber at 150,000, and say where Huber's w lands and why.

Hints:

• 2.2.2's np.where(cond, a, b) builds the charge listing by listing from a mask of the residuals inside the threshold.
• For the two ends, compare the threshold with the largest residual any candidate on the grid leaves, and rewrite delta * (size - delta / 2) as a fixed multiple of size minus a fixed amount.

import numpy as np
import pandas as pd

listings = pd.read_csv("listings.csv")
all_prices = listings["price"].to_numpy()
candidates = np.arange(0, 1_000_001, 1_000)


def huber(residuals, delta):
    size = np.abs(residuals)
    # half the square inside the threshold, the straight line past it
    return np.where(size <= delta, residuals ** 2 / 2, delta * (size - delta / 2))


def scan(charge):
    # the mean charge of every candidate, and the lowest candidate
    scores = np.array([np.mean(charge(all_prices - c)) for c in candidates])
    return candidates[np.argmin(scores)], scores


print(np.median(all_prices), all_prices.mean())
for delta in [1_000, 50_000, 150_000, 500_000, 2_000_000]:
    best, scores = scan(lambda r: huber(r, delta))
    past = np.sum(np.abs(all_prices - best) > delta)
    print(delta, best, past)

337000.0 380763.0
1000 337000 2974
50000 339000 2273
150000 348000 1057
500000 377000 71
2000000 381000 0

# 1930000, the largest residual any candidate leaves
print(np.abs(all_prices - candidates[:, np.newaxis]).max())
best_sq, sq_scores = scan(lambda r: r ** 2)
best_wide, wide_scores = scan(lambda r: huber(r, 2_000_000))
# True: every candidate's Huber score is half its squared score
print(np.allclose(wide_scores, sq_scores / 2))

has_area = listings.dropna(subset=["area_sqm"])
area = has_area["area_sqm"].to_numpy()
price = has_area["price"].to_numpy()
grid = np.arange(0, 6_001, 1.0)
for name, charge in [("squared", lambda r: r ** 2),
                     ("absolute", np.abs),
                     ("huber", lambda r: huber(r, 150_000))]:
    scores = np.array([np.mean(charge(price - w * area)) for w in grid])
    print(name, grid[np.argmin(scores)])

squared 3750.0
absolute 3634.0
huber 3667.0

Additional Resources

Optional, and nothing later depends on either. Worth opening if you want quantiles other than the median, or an estimator built on the Assignment's loss.

• Official docs - scikit-learn, "Quantile regression" example, the "Fitting a QuantileRegressor" section. Three QuantileRegressor lines at the 5%, 50% and 95% quantiles, with the points outside the central 90% band marked, on a target whose spread grows with the feature and on one with a long upper tail.
• Official docs - scikit-learn User Guide, "Linear Models", the "Huber Regression" section. HuberRegressor, a linear model fitted under Huber loss, whose threshold epsilon is measured against a scale the model estimates along with its weights.

Summary

A fit answers the question its loss asks, and the question is easy to miss because it is never printed. Squared error asks for the price that keeps the squared misses smallest, and the answer is the mean. Absolute error asks for the price that keeps the plain misses smallest, and the answer is the median. The model, the data and the call that fits it can be identical, and the loss alone decides which center the parameters chase.

What separates the two is how a listing's distance counts. Under absolute error only its side counts, so a price raised tenfold by a typo leaves the median line exactly where it was. Under squared error each distance is multiplied by itself, so the listings a fit misses most carry a share of the loss far beyond their number. Neither behavior is a flaw, each is right for some question, and the long upper tail of prices is where the choice shows.

The loss also belongs to the fit and not to the report. It is chosen for a search that needs one bottom and a reliable way down to it, which is why LinearRegression runs on squared error, with its single point and its one-step projection. The number a person reads is chosen afterward, for that person's question, in currency and on rows the fit never saw, and keeping the two jobs apart is what turns a printed number in squared currency into a statement a pricing team can act on.

Cheat Sheet

Part A - Reference Table

Part B - Compact Code Block

import numpy as np
from sklearn.linear_model import LinearRegression, QuantileRegressor

grid = np.arange(0, 6_001, 1.0)
losses = np.array([np.mean((y - w * x) ** 2) for w in grid])
best = grid[np.argmin(losses)]  # check it is not an edge
np.mean((y - c) ** 2)           # lowest at c = y.mean()
np.mean(np.abs(y - c))          # lowest at c = np.median(y)
LinearRegression().fit(X, y)    # squared error
# absolute error, the median line
QuantileRegressor(quantile=0.5, alpha=0.0).fit(X, y)
# report the test error in currency
np.abs(y_test - model.predict(X_test)).mean()