Demo: Which Curve Made This Data?

A month at the ice-cream stand: each day we record the temperature \(x\) and how many cones we sold. Nature generated these points from a hidden function plus noise, \[y = f(x) + \eta,\] and nature is not telling us \(f\). All we get is the scatter.

import numpy as np
import matplotlib.pyplot as plt
import scienceplots

plt.style.use(["science", "no-latex"])
TEAL, CARDINAL, GRAY = "#009090", "#9c1b33", "#c9c9c9"
rng = np.random.default_rng(5)

sigma = 5                      # nature's noise level
temps = rng.uniform(58, 100, size=30)
hidden_f = lambda t: 60 - 0.03 * (t - 85) ** 2      # nature's secret
sales = hidden_f(temps) + rng.normal(0, sigma, size=temps.shape)

fig, ax = plt.subplots(figsize=(7, 3))
ax.scatter(temps, sales, s=18, color=GRAY)
ax.set_xlabel("Temperature (°F)")
ax.set_ylabel("Cones sold")
plt.show()

Three theories walk into the shop:

• the intern’s: sales grow linearly with temperature (a straight line);
• mine: sales climb, peak in the mid-80s, then sag when it is too hot to leave the house (a gentle arc);
• a rival stand owner’s: an elaborate 13-parameter curve that threads through the points.

Which curve most plausibly generated this data?

ts = np.linspace(58, 100, 300)
line = lambda t: 0.36 * t + 17
arc = lambda t: 60 - 0.03 * (t - 85) ** 2
wiggle = np.polynomial.Polynomial.fit(temps, sales, 13)

fig, ax = plt.subplots(figsize=(7, 3))
ax.scatter(temps, sales, s=18, color=GRAY)
ax.plot(ts, line(ts), color="black", linewidth=1.4, label="Intern's line")
ax.plot(ts, arc(ts), color=TEAL, linewidth=1.6, label="My arc")
ax.plot(ts, wiggle(ts), color=CARDINAL, linewidth=1.2, label="Rival's wiggle")
ax.set_xlabel("Temperature (°F)")
ax.set_ylabel("Cones sold")
ax.set_ylim(28, 75)
ax.legend(frameon=False)
plt.show()

How do we judge “plausibly”? If a theory is right, then each day’s gap between the curve and the point is pure noise \(\eta\). A theory that requires huge gaps asks us to believe in a wildly unlikely run of noise. So as a first, crude plausibility score, add up the squared gaps; a smaller total is more plausible:

def score(curve, t, y):
    return np.mean((curve(t) - y) ** 2)

print("mean squared gap, this month:")
print(f"  intern's line   {score(line, temps, sales):8.1f}")
print(f"  my arc          {score(arc, temps, sales):8.1f}")
print(f"  rival's wiggle  {score(wiggle, temps, sales):8.1f}")

# a fresh month of days from the same hidden process
temps2 = rng.uniform(58, 100, size=30)
sales2 = hidden_f(temps2) + rng.normal(0, sigma, size=temps2.shape)

print("\nmean squared gap, NEXT month:")
print(f"  intern's line   {score(line, temps2, sales2):8.1f}")
print(f"  my arc          {score(arc, temps2, sales2):8.1f}")
print(f"  rival's wiggle  {score(wiggle, temps2, sales2):8.1f}")

mean squared gap, this month:
  intern's line       90.7
  my arc              26.8
  rival's wiggle       9.8

mean squared gap, NEXT month:
  intern's line      117.6
  my arc              19.2
  rival's wiggle     220.2

On this month’s data the rival’s wiggle wins: it threads the points, so its gaps are smallest. But on a fresh month from the same process the wiggle falls apart while the arc’s score barely moves. The wiggle fit the noise in the month it was shown, not the process behind it. (A whole lecture, Methodology, is devoted to this trap.)

And here is the reveal: nature’s hidden function was the arc, \(f(x) = 60 - 0.03(x - 85)^2\). Yet even the true curve misses every day by about \(\sigma = 5\) cones; its mean squared gap hovers near \(\sigma^2 = 25\), and no curve can beat that on fresh data.

What if nature were noisier?

Turn the one knob nature has: \(\sigma\). Try sigma = 0 (a world with no birthday parties) and sigma = 15 (chaos).

sigma_new = 15    # change me!

temps3 = rng.uniform(58, 100, size=30)
sales3 = hidden_f(temps3) + rng.normal(0, sigma_new, size=temps3.shape)

fig, ax = plt.subplots(figsize=(7, 3))
ax.scatter(temps3, sales3, s=18, color=GRAY)
ax.plot(ts, hidden_f(ts), color=TEAL, linewidth=1.6, label="The true curve $f$")
ax.set_xlabel("Temperature (°F)")
ax.set_ylabel("Cones sold")
ax.legend(frameon=False)
plt.show()

print(f"true curve's mean squared gap: {score(hidden_f, temps3, sales3):.1f}"
      f"   (sigma^2 = {sigma_new**2})")

true curve's mean squared gap: 227.9   (sigma^2 = 225)

Punchline: even the true function misses by the noise, so zero error is not the goal, and a model that reaches it has fit the day-to-day luck instead of the pattern. The best any predictor can do is the conditional mean \(\mathbb{E}[Y \mid X = x]\), and Problem 5 proves its error floor \(\mathbb{E}[\mathrm{Var}(Y \mid X)]\) exactly. Next lecture, today’s crude “sum of squared gaps” score stops being crude: it falls out of maximum likelihood as a theorem.