Museum regression results

A read-only report on the regression model that the museum_analytics FastAPI server has loaded. Everything here comes from the API over HTTP. The notebook never calls /train, so it shows the model the server built from its own data at startup and leaves it unchanged.

Endpoint	Used for
`GET /health`	check the server is up and has a model
`GET /model`	fit statistics and the fitted equation
`GET /data`	every museum the model was fitted on, with its prediction
`POST /predict`	the fitted curve, and predictions for new cities
`GET /residuals`	museums that most beat or miss what their city size predicts

In [2]:

import os
import requests
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from matplotlib.ticker import FuncFormatter, LogLocator, NullFormatter

# Inside docker-compose this is http://museum_analytics_app:8000; outside it, localhost.
API = os.getenv("MUSEUM_ANALYTICS_URL", "http://localhost:8000")


def get(path, **params):
    r = requests.get(f"{API}{path}", params=params, timeout=10)
    r.raise_for_status()
    return r.json()


def post(path, payload):
    r = requests.post(f"{API}{path}", json=payload, timeout=10)
    r.raise_for_status()
    return r.json()


BLUE, ORANGE, INK, MUTED, GRID = "#2a78d6", "#eb6834", "#0b0b0b", "#52514e", "#e6e5e0"
plt.rcParams.update({
    "axes.spines.top": False, "axes.spines.right": False,
    "axes.edgecolor": MUTED, "axes.labelcolor": INK, "axes.titlecolor": INK,
    "xtick.color": MUTED, "ytick.color": MUTED,
    "axes.grid": True, "axes.axisbelow": True, "grid.color": GRID, "grid.linewidth": 0.8,
    "legend.frameon": False, "figure.dpi": 110,
})


def human(x, _=None):
    """1_500_000 -> '1.5M'"""
    for div, suf in ((1e9, "B"), (1e6, "M"), (1e3, "K")):
        if abs(x) >= div:
            return f"{x / div:.3g}{suf}"
    return f"{x:.3g}"


def human_axes(ax, log=False):
    """Readable tick labels; on log axes tick at 1, 2, 5 x 10^k."""
    for axis in (ax.xaxis, ax.yaxis):
        if log:
            axis.set_major_locator(LogLocator(base=10, subs=(1, 2, 5)))
            axis.set_minor_formatter(NullFormatter())
        axis.set_major_formatter(FuncFormatter(human))

1. Check the server

In [3]:

try:
    health = get("/health")
except requests.ConnectionError:
    raise SystemExit(f"Can't reach the API at {API}. Start it with: docker compose up --build -d")

if not health["model_loaded"]:
    raise SystemExit("The server is up but has no model loaded. POST /train first.")
health

Out [3]:

{'status': 'ok', 'model_loaded': True}

2. Model summary

GET /model returns the fit statistics. With the default log scale the model is log10(visitors) ~ log10(population), so the slope is an elasticity: how visitors scale with city size.

In [4]:

summary = get("/model")
is_log = summary["scale"] == "log"

stats = pd.DataFrame({
    "value": [summary["n_samples"], summary["scale"], summary["aggregate_by_city"],
              summary["slope"], summary["intercept"], summary["r2"],
              summary["cv_r2_5fold"], summary["pearson_r"], summary["spearman_rho"]],
    "meaning": ["rows fitted", "axis transform", "one row per city?",
                "change in (log) visitors per unit of (log) population", "fitted intercept",
                "share of variance explained (training data)",
                "R² on held-out folds (5-fold CV)",
                "linear correlation", "rank correlation"],
}, index=["n", "scale", "aggregate_by_city", "slope", "intercept", "R²",
          "CV R²", "Pearson r", "Spearman ρ"])

print(summary["equation"])
print(summary["interpretation"])
fit = summary["cv_r2_5fold"] if summary["cv_r2_5fold"] is not None else summary["r2"]
if fit < 0.3:
    print(f"\nCaution: city population explains only about {max(fit, 0):.0%} of the variation in "
          "visitors, so treat predictions as rough baselines.")
stats

Out [4]:

visitors = 338,577.3 * population^0.128
A 10x larger city is associated with 1.34x the visitors.

Caution: city population explains only about 2% of the variation in visitors, so treat predictions as rough baselines.

	value	meaning
n	71	rows fitted
scale	log	axis transform
aggregate_by_city	False	one row per city?
slope	0.128263	change in (log) visitors per unit of (log) pop...
intercept	5.529658	fitted intercept
R²	0.084989	share of variance explained (training data)
CV R²	0.016122	R² on held-out folds (5-fold CV)
Pearson r	0.291528	linear correlation
Spearman ρ	0.331742	rank correlation

3. The fitted line

Each dot is a museum. The line is drawn by sending a grid of populations to POST /predict, so it is exactly what the server would answer. The five museums furthest from the line are labelled.

In [5]:

points = pd.DataFrame(get("/data"))

grid = (np.geomspace if is_log else np.linspace)(points.Population.min(), points.Population.max(), 200)
curve = pd.DataFrame(post("/predict", {"populations": grid.tolist()}))

# Residual on the scale the model was fitted on
if is_log:
    points["residual"] = np.log10(points.Visitors / points.Predicted)
else:
    points["residual"] = points.Visitors - points.Predicted

fig, ax = plt.subplots(figsize=(9, 6))
ax.scatter(points.Population, points.Visitors, s=40, color=BLUE, alpha=0.75,
           edgecolor="white", linewidth=1, label="Museum", zorder=3)
ax.plot(curve.population, curve.predicted_visitors, color=ORANGE, lw=2,
        label=f"Fit: {summary['equation']}", zorder=4)

for _, r in points.loc[points.residual.abs().nlargest(5).index].iterrows():
    ax.annotate(r.Name, (r.Population, r.Visitors), xytext=(6, 4),
                textcoords="offset points", fontsize=8, color=MUTED)

if is_log:
    ax.set_xscale("log")
    ax.set_yscale("log")
human_axes(ax, is_log)
ax.set_xlabel("City population")
ax.set_ylabel("Annual visitors")
cv = summary["cv_r2_5fold"]
ax.set_title(f"Visitors vs. city population    R² = {summary['r2']:.2f}"
             + (f"  (cross-validated {cv:.2f})" if cv is not None else ""),
             loc="left", fontsize=12)
ax.legend(loc="upper left")
plt.tight_layout()
plt.show()

4. How good are the predictions?

Left: actual vs. predicted visitors. Points on the dashed line are predicted perfectly; above it the museum beats its prediction. Right: the spread of residuals. A narrow, centred spread means city size explains a lot; a wide one means other factors matter more.

In [6]:

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 5))

ax1.scatter(points.Predicted, points.Visitors, s=40, color=BLUE, alpha=0.75,
            edgecolor="white", linewidth=1, zorder=3)
lo = min(points.Predicted.min(), points.Visitors.min()) * 0.9
hi = max(points.Predicted.max(), points.Visitors.max()) * 1.1
ax1.plot([lo, hi], [lo, hi], color=MUTED, lw=1, ls="--", label="Perfect prediction")
if is_log:
    ax1.set_xscale("log")
    ax1.set_yscale("log")
ax1.set_xlim(lo, hi)   # same range on both axes so the diagonal is a true 45° line
ax1.set_ylim(lo, hi)
ax1.set_aspect("equal")
human_axes(ax1, is_log)
ax1.set_xlabel("Predicted visitors")
ax1.set_ylabel("Actual visitors")
ax1.set_title("Actual vs. predicted", loc="left", fontsize=12)
ax1.legend(loc="lower right")

ax2.hist(points.residual, bins=20, color=BLUE, edgecolor="white", linewidth=2)
ax2.axvline(0, color=MUTED, lw=1)
if is_log:
    ticks = np.array([0.1, 0.2, 0.5, 1, 2, 5, 10])
    ticks = ticks[(np.log10(ticks) >= points.residual.min() - 0.1) &
                  (np.log10(ticks) <= points.residual.max() + 0.1)]
    ax2.set_xticks(np.log10(ticks), [f"{t:g}×" for t in ticks])
    ax2.set_xlabel("Actual ÷ predicted visitors (log scale)")
else:
    ax2.xaxis.set_major_formatter(FuncFormatter(human))
    ax2.set_xlabel("Actual − predicted visitors")
ax2.set_ylabel("Museums")
ax2.set_title("Residual distribution", loc="left", fontsize=12)
ax2.grid(axis="x", visible=False)

plt.tight_layout()
plt.show()

within_2x = (points.residual.abs() <= np.log10(2)).mean() if is_log else None
if within_2x is not None:
    print(f"{within_2x:.0%} of museums are within 2× of their predicted visitors.")

82% of museums are within 2× of their predicted visitors.

5. Which museums beat (or miss) their city's size?

GET /residuals ranks museums by actual ÷ predicted visitors. Above 1× (blue) the museum draws more than its city's size alone would suggest; below 1× (orange) it draws fewer.

In [7]:

N = 8
over = pd.DataFrame(get("/residuals", top=N, order="over"))
under = pd.DataFrame(get("/residuals", top=N, order="under"))
res = (pd.concat([over, under.iloc[::-1]])
       .drop_duplicates("Name")
       .sort_values("Ratio"))

fig, ax = plt.subplots(figsize=(9, 0.35 * len(res) + 1.2))
colors = [BLUE if r >= 1 else ORANGE for r in res.Ratio]
ax.barh(res.Name + " (" + res.City + ")", np.log10(res.Ratio), color=colors, height=0.7)
ax.axvline(0, color=MUTED, lw=1)
ticks = np.array([0.05, 0.1, 0.2, 0.5, 1, 2, 5, 10, 20])
ticks = ticks[(np.log10(ticks) >= np.log10(res.Ratio).min() - 0.1) &
              (np.log10(ticks) <= np.log10(res.Ratio).max() + 0.1)]
ax.set_xticks(np.log10(ticks), [f"{t:g}×" for t in ticks])
ax.set_xlabel("Actual ÷ predicted visitors (log scale)")
ax.set_title(f"Top {N} over- and under-performers", loc="left", fontsize=12)
ax.grid(axis="y", visible=False)
plt.tight_layout()
plt.show()

res[["Name", "City", "Country", "Population", "Visitors", "Predicted", "Ratio"]] \
    .sort_values("Ratio", ascending=False).round(2).reset_index(drop=True)

Out [7]:

	Name	City	Country	Population	Visitors	Predicted	Ratio
0	Louvre	Paris	France	11060000.0	9000000	2710895.23	3.32
1	Shenzhen Museum	Shenzhen	China	14678000.0	6805000	2811109.70	2.42
2	National Museum of China	Beijing	China	21893095.0	7031700	2959030.62	2.38
3	British Museum	London	United Kingdom	11262000.0	6440120	2717195.79	2.37
4	Natural History Museum, South Kensington	London	United Kingdom	11262000.0	6301972	2717195.79	2.32
5	China Science and Technology Museum	Beijing	China	21893095.0	6421000	2959030.62	2.17
6	Nanjing Museum	Nanjing	China	9341685.0	5680000	2652816.45	2.14
7	Metropolitan Museum of Art	New York City	United States	19268388.0	5984091	2910956.74	2.06
8	Museum of Science	Boston	United States	4453352.0	1324000	2412348.33	0.55
9	Palacio de Cristal del Retiro	Madrid	Spain	6211000.0	1318823	2517507.58	0.52
10	Smithsonian American Art Museum (with Renwick ...	Washington, D.C.	United States	5230370.0	1273450	2462626.87	0.52
11	Art Institute of Chicago	Chicago	United States	8609571.0	1324241	2625192.03	0.50
12	Getty Center	Los Angeles	United States	11984083.0	1301332	2738940.93	0.48
13	Centro Cultural Banco do Brasil	São Paulo	Brazil	23086000.0	1364208	2979235.54	0.46
14	Chinese Aviation Museum	Beijing	China	21893095.0	1292278	2959030.62	0.44
15	Moscow Kremlin Museum	Moscow	Russia	19100000.0	1240113	2907681.34	0.43

6. Predict for new cities

Edit cities and re-run. A single city uses GET /predict?population=...; a batch uses POST /predict. The predictions are drawn on the fitted line from section 3.

In [8]:

cities = {"Reykjavík": 140_000, "Lyon": 520_000, "Montréal": 1_780_000,
          "Mexico City": 9_200_000, "Tokyo": 14_000_000}

single = get("/predict", population=cities["Montréal"])
print(f"Montréal (single call): {human(single['predicted_visitors'])} visitors")

pred = pd.DataFrame(post("/predict", {"populations": list(cities.values())}), index=list(cities))

fig, ax = plt.subplots(figsize=(9, 6))
ax.scatter(points.Population, points.Visitors, s=30, color=GRID, edgecolor=MUTED,
           linewidth=0.5, label="Fitted museums", zorder=2)
ax.plot(curve.population, curve.predicted_visitors, color=ORANGE, lw=2, label="Fit", zorder=3)
ax.scatter(pred.population, pred.predicted_visitors, s=70, color=BLUE,
           edgecolor="white", linewidth=2, label="Predicted city", zorder=4)
for city, r in pred.iterrows():
    ax.annotate(f"{city}: {human(r.predicted_visitors)}", (r.population, r.predicted_visitors),
                xytext=(8, -12), textcoords="offset points", fontsize=9, color=INK)
if is_log:
    ax.set_xscale("log")
    ax.set_yscale("log")
human_axes(ax, is_log)
ax.set_xlabel("City population")
ax.set_ylabel("Predicted annual visitors per museum")
ax.set_title("Predicted visitors for new cities", loc="left", fontsize=12)
ax.legend(loc="upper left")
plt.tight_layout()
plt.show()

pred.rename(columns={"population": "Population", "predicted_visitors": "Predicted visitors"}).round(0)

Out [8]:

Montréal (single call): 2.14M visitors

	Population	Predicted visitors
Reykjavík	140000.0	1547807.0
Lyon	520000.0	1831517.0
Montréal	1780000.0	2144653.0
Mexico City	9200000.0	2647621.0
Tokyo	14000000.0	2794109.0

347 KiB Raw Blame History