kuh-no-bowls/knb_stats/writeup.ipynb

Reading a Theme Park's Mind with PCA

What 600,000 queue-time readings can tell you about how people move through Knoebels

I'm working my way through An Introduction to Statistical Learning and Gilbert Strang's linear algebra lectures, and I got impatient. Principal Component Analysis sits a few chapters ahead of where I am, but I had a dataset I couldn't stop thinking about, so I skipped ahead and used it as an excuse to learn PCA by doing. This notebook is me figuring it out live.

The data. Every few minutes, a scraper records the posted wait time for each ride at Knoebels (a wonderful, weird old amusement park in Pennsylvania). Two seasons of that adds up to ~600k readings across ~57 rides.

The question. If I don't tell the algorithm anything about the rides — not which are roller coasters, not which are for toddlers, not where they sit in the park — can it still discover the hidden "axes" that explain how the whole park behaves? That's exactly what PCA is for: it finds the handful of directions that explain the most variation in the data, and leaves us to interpret them.

Spoiler: it finds three, and they turn out to be how busy the park is, which half of the park you're in, and what time of day it is — none of which it was ever told.

The data: one big table of waits

The natural shape for this is a table where every row is a moment in time (a 5-minute snapshot of the park) and every column is a ride. Each cell is the wait, in minutes, for that ride at that moment.

A couple of small judgment calls in building it:

• Forward-fill, then zero. If a ride didn't report at 2:05 we carry its 2:00 value forward (the line probably didn't teleport). Anything still missing — like early morning before the ride opens — becomes 0.
• Drop powersurge. It hasn't been running for most of the time since I began collecting data; ~90% of its column is imputed zeros, so it only adds noise.

In [24]:

import os
from pathlib import Path

import pandas as pd
from dotenv import load_dotenv
from sqlalchemy import create_engine, text

load_dotenv(Path.cwd() / ".env")
engine = create_engine(
    "postgresql+psycopg2://{user}:{pw}@{host}:{port}/{db}".format(
        user=os.environ["PG_USER"], pw=os.environ["PG_PASSWORD"],
        host=os.environ["PG_HOST"], port=os.environ["PG_PORT"], db=os.environ["PG_DB"],
    )
)

df = pd.read_sql(
    text(
        """
        SELECT a.time_stamp,
               r.ride        AS ride_name,
               a.queue_time  AS queue_time_sec
          FROM knoebels."LZ_attractions_io" a
          JOIN knoebels.ride_master r ON r.attractions_dot_io_id = a._id
         WHERE a.queue_time IS NOT NULL
        """
    ),
    engine,
)
df["queue_time_min"] = df["queue_time_sec"] / 60.0
df["time_stamp"] = pd.to_datetime(df["time_stamp"])


def ride_time_matrix(df, freq="5min"):
    """Long table -> (time x ride) matrix of mean wait in minutes per time bin."""
    snapped = df.assign(ts=df["time_stamp"].dt.floor(freq))
    return (
        snapped.pivot_table(index="ts", columns="ride_name",
                            values="queue_time_min", aggfunc="mean")
        .sort_index()
        .ffill()      # carry the last known wait forward
        .fillna(0)    # anything still missing (pre-open) = no line
    )


df_matrix = ride_time_matrix(df).drop(columns=["powersurge"], errors="ignore")
print(f"{len(df):,} readings  ->  matrix of {df_matrix.shape[0]:,} time-snapshots x {df_matrix.shape[1]} rides")
#df_matrix.iloc[:3, :5]
df_matrix.info()

657,498 readings  ->  matrix of 12,535 time-snapshots x 56 rides
<class 'pandas.DataFrame'>
DatetimeIndex: 12535 entries, 2024-05-15 12:20:00 to 2026-05-31 17:30:00
Data columns (total 56 columns):
 #   Column              Non-Null Count  Dtype  
---  ------              --------------  -----  
 0   antique_cars        12535 non-null  float64
 1   balloon_race        12535 non-null  float64
 2   black_diamond       12535 non-null  float64
 3   bumper_cars         12535 non-null  float64
 4   cosmotron           12535 non-null  float64
 5   downdraft           12535 non-null  float64
 6   fandango            12535 non-null  float64
 7   flyer               12535 non-null  float64
 8   flying_tigers       12535 non-null  float64
 9   flying_turns        12535 non-null  float64
 10  galleon             12535 non-null  float64
 11  giant_flume         12535 non-null  float64
 12  giant_wheel         12535 non-null  float64
 13  goin'_buggy         12535 non-null  float64
 14  grand_carousel      12535 non-null  float64
 15  hand_cars           12535 non-null  float64
 16  haunted_mansion     12535 non-null  float64
 17  helicopters         12535 non-null  float64
 18  impulse             12535 non-null  float64
 19  italian_trapeze     12535 non-null  float64
 20  jet_skyfighter      12535 non-null  float64
 21  kiddie_boats        12535 non-null  float64
 22  kiddie_bumper_cars  12535 non-null  float64
 23  kiddie_firetrucks   12535 non-null  float64
 24  kiddie_himalaya     12535 non-null  float64
 25  kiddie_wheel        12535 non-null  float64
 26  kiddie_whip         12535 non-null  float64
 27  kozmo's_kurves      12535 non-null  float64
 28  looper              12535 non-null  float64
 29  merry_mixer         12535 non-null  float64
 30  motor_boats         12535 non-null  float64
 31  ole_smokey          12535 non-null  float64
 32  panther_cars        12535 non-null  float64
 33  paradrop            12535 non-null  float64
 34  paratrooper         12535 non-null  float64
 35  pete's_fleet        12535 non-null  float64
 36  phoenix             12535 non-null  float64
 37  pioneer_train       12535 non-null  float64
 38  pony_carts          12535 non-null  float64
 39  red_baron           12535 non-null  float64
 40  ribbit              12535 non-null  float64
 41  roto_jets           12535 non-null  float64
 42  s&g_carousel        12535 non-null  float64
 43  satellite           12535 non-null  float64
 44  scenic_skyway       12535 non-null  float64
 45  sklooosh            12535 non-null  float64
 46  spanish_bambini     12535 non-null  float64
 47  stratosfear         12535 non-null  float64
 48  super_round-up      12535 non-null  float64
 49  tea_cups            12535 non-null  float64
 50  tilt-a-whirl        12535 non-null  float64
 51  tornado             12535 non-null  float64
 52  tumbling_timbers    12535 non-null  float64
 53  twister             12535 non-null  float64
 54  umbrella_ride       12535 non-null  float64
 55  whipper             12535 non-null  float64
dtypes: float64(56)
memory usage: 5.5 MB

Step 1: put every ride on the same footing

Here's the first thing that bit me. PCA hunts for directions of maximum variance — but "variance" is measured in whatever units the data is in. A headline roller coaster might swing between 5 and 45 minutes; a kiddie ride swings between 0 and 6. If I feed in raw minutes, PCA will conclude the coaster is "more important" purely because its numbers are bigger, and the small rides will be invisible.

The fix is to standardize every column: subtract its mean and divide by its standard deviation, so each ride has mean 0 and variance 1. Now a "big day" for the kiddie carousel counts just as much as a big day for the coaster. We're asking about patterns of co-movement, not raw magnitudes.

In [17]:

from sklearn.preprocessing import StandardScaler

scaler = StandardScaler()
scaled = scaler.fit_transform(df_matrix)   # each column now mean 0, variance 1
pd.DataFrame(scaled, columns=df_matrix.columns).describe().loc[["mean", "std"]].iloc[:, :4].round(2)

Out [17]:

ride_name	antique_cars	balloon_race	black_diamond	bumper_cars
mean	0.0	-0.0	-0.0	-0.0
std	1.0	1.0	1.0	1.0

Step 2: ask for the axes of variation

This is the part I'd been intimidated by, and it turned out to be a clean idea.

We have ~57 ride-columns, but they're hugely redundant — when the park is mobbed, everything gets busy together. PCA asks: can I describe most of that 57-dimensional wobble with just a few new axes? Each axis (a principal component) is a weighted blend of the original rides, chosen so the first one captures as much variance as possible, the second captures as much of what's left as possible while being perpendicular to the first, and so on.

The linear algebra underneath is the singular value decomposition. Our standardized matrix X factors as

X = U\,\Sigma\,V^{\top}.

The columns of V are the new axes — the loadings, telling us how much each ride contributes to each component. U\Sigma gives the scores — where each time-snapshot lands along those axes. And the singular values in \Sigma, squared, tell us how much variance each axis explains. sklearn's PCA does exactly this and hands back the pieces with friendly names.

I asked for the top 3 components. Look at how concentrated the signal is:

In [18]:

from sklearn.decomposition import PCA

pca = PCA(n_components=3)
principal_components = pca.fit_transform(scaled)

scores = pd.DataFrame(principal_components, index=df_matrix.index,
                      columns=["PC1", "PC2", "PC3"])
loadings = pd.DataFrame(pca.components_.T, index=df_matrix.columns,
                        columns=["PC1", "PC2", "PC3"])

for i, pc in enumerate(["PC1", "PC2", "PC3"]):
    print(f"{pc}: {pca.explained_variance_ratio_[i]:5.1%} of all variance")
print(f"{'':4} {pca.explained_variance_ratio_[:3].sum():5.1%} captured by just 3 axes (out of 56)")

PC1: 40.4% of all variance
PC2: 13.0% of all variance
PC3:  6.7% of all variance
     60.1% captured by just 3 axes (out of 56)

Three axes out of fifty-six explain ~60% of everything. The rest is ride-specific noise. Now the fun part — what are these three axes? PCA won't tell us; it just hands us numbers. We have to read them.

My approach for each component: look at which rides load heavily on it (the loadings), form a hypothesis about what real-world thing that represents, then validate the hypothesis against something the PCA never saw.

PC1 — the whole park breathes together

The first component is the easiest. Every single ride loads on it with the same sign — they all move together. That's the signature of one global force: overall busyness. Empty Tuesday morning, the whole park is near zero; packed Saturday afternoon, everything spikes at once.

The test: if PC1 is really "busyness," then a snapshot's PC1 score should track the plain average wait across all rides at that moment. To make that easy to eyeball, I put both quantities on the same standardized scale — so if PC1 were exactly the average wait, every dot would land on the dashed 45° line. The red line is the actual relationship (its slope is the correlation). For PC1 the red line sits almost on top of the diagonal:

In [19]:

import numpy as np
import matplotlib.pyplot as plt


def validation_plot(real_quantity, pc_score, score_name, xlabel):
    """Put a PC score and the real-world quantity it should track on the same
    standardized scale. If the PC *were* that quantity every dot would sit on the
    dashed 45° line; the red line is the actual fit (its slope is the correlation
    r). The closer red hugs the diagonal, the stronger the validation."""
    zx = (real_quantity - real_quantity.mean()) / real_quantity.std()
    zy = (pc_score - pc_score.mean()) / pc_score.std()
    r = real_quantity.corr(pc_score)

    fig, ax = plt.subplots(figsize=(7, 6.5))
    ax.scatter(zx, zy, s=4, alpha=0.12, color="steelblue")
    line = np.array([-3, 3])
    ax.plot(line, line, "k--", lw=1.3, label="perfect match (y = x)")
    ax.plot(line, r * line, color="crimson", lw=2.5,
            label=f"actual relationship (r = {r:.2f})")
    ax.set(xlim=(-3, 3), ylim=(-3, 3), aspect="equal",
           xlabel=f"{xlabel}  (standardized)",
           ylabel=f"{score_name} score (standardized)")
    ax.set_title(f"{score_name}: does it track the real-world quantity?")
    ax.legend(loc="upper left")
    plt.tight_layout()


park_mean = df_matrix.mean(axis=1)
validation_plot(park_mean, scores["PC1"], "PC1", "Park-wide mean wait")

PC2 — Kiddieland vs. Everything else

PC2 is a contrast: some rides load positive, others negative. Reading the extremes, the negative end is dominated by the tiny-tots rides (kiddie_boats, pony_carts, red_baron...) and the positive end by the big rides in the rest of the park. So my first guess was "kiddie rides vs. everything else."

The category tag is a trap

The park's own data tags each ride Kiddie / Family / Thrill, so I reached for that to define the two groups — and got a muddy result. Two problems:

1. The tag values secretly carry a trailing space ("Kiddie ", not "Kiddie"), so naive equality checks silently matched nothing. (Lesson learned: a value_counts() printout will never show you that.)
2. Even cleaned up, the tag misfiles rides relative to how they behave. The big Ferris wheel (giant_wheel) and the large grand_carousel are tagged Family but behave like main-midway rides; the small s&g_carousel is tagged Family but behaves like a Kiddieland ride.

PC2 isn't tracking the taxonomy — it's tracking physical geography: the Kiddieland corner vs. the main midway. So let me define the two groups from geometry the PCA never touched, and see if the contrast gets sharper.

Validating from the map

Each ride has a GPS coordinate. I run k-means on just the lat/lon (no queue data, no PCA) to carve the park into areas, and call the cluster with the lowest mean PC2 loading "Kiddieland." Then I build the contrast main-midway mean − Kiddieland mean and check it against PC2.

(Small but important detail: standardize the coordinates first. On raw degrees the clustering is unstable from one random seed to the next; standardized, it's identical every run.)

In [20]:

from sklearn.cluster import KMeans

# category tag, with the sneaky trailing space trimmed off
category = pd.read_sql(
    text(
        """
        SELECT r.ride AS ride_name, trim(c.name) AS category
          FROM knoebels.ride_master r
          JOIN knoebels."LZ_attractions_io_poi" p ON p._id = r.attractions_dot_io_id
          LEFT JOIN knoebels."LZ_attractions_io_categories" c ON c._id = p.category
        """
    ),
    engine,
).set_index("ride_name")["category"]

coords = pd.read_sql(
    text(
        """
        SELECT r.ride AS ride_name,
               (p.location)[0]::float AS lat, (p.location)[1]::float AS lon
          FROM knoebels.ride_master r
          JOIN knoebels."LZ_attractions_io_poi" p ON p._id = r.attractions_dot_io_id
         WHERE p.location IS NOT NULL
        """
    ),
    engine,
).set_index("ride_name")

annotated = loadings.join(category.rename("category"))
located = annotated.join(coords, how="inner")

# k-means on STANDARDIZED coordinates only -> reproducible spatial areas
located["area"] = KMeans(n_clusters=5, n_init=10, random_state=0).fit_predict(
    StandardScaler().fit_transform(located[["lat", "lon"]])
)
kiddie_area = located.groupby("area")["PC2"].mean().idxmin()
kiddie_rides = located.index[located["area"] == kiddie_area]
main_rides   = located.index[located["area"] != kiddie_area]

# Does the map agree? Left: each ride at its real GPS point, colored by PC2
# loading. Right: areas from clustering lat/lon alone — no queue data, no PCA.
fig, axes = plt.subplots(1, 2, figsize=(15, 7))
sc = axes[0].scatter(located["lon"], located["lat"], c=located["PC2"],
                     cmap="RdBu_r", s=140, edgecolor="black", lw=0.4)
axes[0].set(title="Each ride at its real location, colored by PC2 loading\n"
                  "(blue = Kiddieland end, red = main-midway end)",
            xlabel="Longitude", ylabel="Latitude")
plt.colorbar(sc, ax=axes[0], label="PC2 loading")
for ride, row in located.iterrows():
    axes[0].annotate(ride, (row["lon"], row["lat"]), fontsize=6, alpha=0.6,
                     xytext=(3, 3), textcoords="offset points")
for area, g in located.groupby("area"):
    label = f"Area {area}" + (" (Kiddieland)" if area == kiddie_area else "")
    axes[1].scatter(g["lon"], g["lat"], label=label, s=140, alpha=0.85)
axes[1].set(title="k-means on standardized lat/lon alone\n(no queue data, no PCA)",
            xlabel="Longitude", ylabel="Latitude")
axes[1].legend()
plt.tight_layout()
print(f"Kiddieland ({len(kiddie_rides)} rides): {', '.join(sorted(kiddie_rides))}")

Kiddieland (19 rides): goin'_buggy, hand_cars, haunted_mansion, helicopters, jet_skyfighter, kiddie_boats, kiddie_bumper_cars, kiddie_firetrucks, kiddie_wheel, kiddie_whip, ole_smokey, panther_cars, pete's_fleet, pony_carts, red_baron, ribbit, s&g_carousel, spanish_bambini, umbrella_ride

The blue (negative-PC2) rides really do huddle in one corner of the park — and clustering the coordinates on their own rediscovers that same corner without ever touching a wait time. So PC2 is genuinely spatial. Now the quantitative check: take main-midway mean − Kiddieland mean at each moment and see how tightly it tracks the PC2 score.

In [21]:

contrast = df_matrix[main_rides].mean(axis=1) - df_matrix[kiddie_rides].mean(axis=1)
validation_plot(contrast, scores["PC2"], "PC2", "Main-midway − Kiddieland wait")

r ≈ 0.78, built from coordinates the model never saw, and giant_wheel / grand_carousel correctly land on the main-midway side while s&g_carousel lands in Kiddieland — exactly the rides the category tag got wrong. PC2 is where in the park you are.

PC3 — morning school groups vs. the evening crowd

PC3 was the one I got wrong at first, so this is the honest version.

My initial hypothesis was ride throughput: "people-eater" walk-on rides vs. low-capacity rides that build a line. It's a tidy story, and it's wrong here:

• The capacity and duration columns I'd have needed are completely empty in this dataset, so I couldn't even compute throughput.
• When I split rides by their PC3 sign and compared actual wait behavior, the two groups were a near-tie (mean wait 9.8 vs. 8.8 min). PC3 simply doesn't separate "builds a line" from "walk-on."

What it does separate is time of day. The PC3 score runs strongly positive at midday and strongly negative after dark. So the rides loading positive are busy in the daytime, and the ones loading negative are busy in the evening — and that maps onto a real thing about how Knoebels operates:

Knoebels runs discounted group rates for schools and daycares, who arrive at open and leave by mid-afternoon. Those big groups of little kids hammer the daytime family flat-rides (tea cups, tilt-a-whirl, paratrooper). After they clear out, the evening belongs to teens and date-night couples on the coasters, bumper cars, and the Ferris wheel at dusk.

The test: de-trend each group's wait by the park-wide average for that hour (so we're looking at shape, not PC1's overall nightly ramp), and watch the two groups pull apart across the day.

In [22]:

daytime_rides = loadings.index[loadings["PC3"] > 0.1]
evening_rides = loadings.index[loadings["PC3"] < -0.1]

hour = df_matrix.index.hour
park_by_hour = df_matrix.mean(axis=1).groupby(hour).mean()

fig, axes = plt.subplots(1, 2, figsize=(14, 5))

# left: the axis itself, by hour
scores.groupby(scores.index.hour)["PC3"].mean().plot(ax=axes[0], marker="o")
axes[0].axhline(0, color="grey", lw=0.5)
axes[0].set(title="PC3 score by hour  (+ daytime / − evening)",
            xlabel="Hour (local)", ylabel="Mean PC3 score")

# right: each group's wait as a share of the park-wide level, de-trended
for name, grp, color in [("daytime rides (PC3 +)", daytime_rides, "tab:orange"),
                         ("evening rides (PC3 −)", evening_rides, "tab:blue")]:
    (df_matrix[grp].mean(axis=1).groupby(hour).mean() / park_by_hour).plot(
        ax=axes[1], marker="o", label=name, color=color)
axes[1].axhline(1.0, color="grey", lw=0.5)
axes[1].set(title="Group wait ÷ park-wide wait, by hour",
            xlabel="Hour (local)", ylabel="Share of park-wide mean wait")
axes[1].legend()
plt.tight_layout()

The daytime group fades through the afternoon while the evening group climbs — and the school-group operating detail is part of the mechanism. PC3 is what time of day it is.

(There's a mild bit of circularity worth flagging: I defined the two groups using the PC3 loadings, so this plot confirms what the axis means more than it's a fully independent test. The genuinely external evidence is the operational fact about group rates.)

What PC3 adds that PC2 doesn't

It's fair to ask what's left for PC3 to do. PC2's score also drifts with the clock — kiddie-heavy in the morning, midway-heavy by night — so a third time-flavored axis can feel redundant. It isn't, and seeing why gets at the heart of why PCA stacks several axes.

The mechanics: each component assigns every ride a single number (its loading), and the axes are built to be independent — here PC2 and PC3 correlate 0.00. PC2 spends its one-number-per-ride on a single question — kiddie corner or main midway? — so every grown-up ride gets a similar positive loading. That makes PC2 structurally blind to the fact that, say, the tea cups and the giant Ferris wheel behave nothing alike. Both are just "midway" to it.

In [ ]:

# PC2 calls these two rides almost identical (both lean "main midway").
# PC3 sends them to opposite ends.
loadings.loc[["tilt-a-whirl", "giant_wheel"], ["PC2", "PC3"]].round(3)

PC2 puts them practically side by side (+0.16 and +0.20); PC3 sends them to opposite poles (+0.15 vs −0.18). And that matches the park exactly: the tilt-a-whirl is a daytime ride the school-and-daycare crowd swarms at opening, while the giant wheel is an evening ride — dusk views, the date-night crowd. To describe a ride as "midway and daytime" versus "midway and evening" takes two numbers, and the second one is PC3.

It even shows up in the shape of each axis across a day: PC2 is a steady all-day drift, PC3 a midday hump. So the three components stack into one clean read on the park:

• PC1 — how many people are here? (overall busyness)
• PC2 — which half of the park are they in? (kiddie corner ↔ main midway)
• PC3 — among the midway rides, which kind is peaking? (daytime family flat-rides ↔ evening marquee rides)

The gut-check for what "independent axes" really buys you: two moments can share an identical PC1 and PC2 — same crowd size, same kiddie/midway balance — and still split on PC3. A midday with the flat rides surging versus an evening with the coasters and bumper cars going: same volume, same spatial tilt, different ride mix. PC3 is the only one of the three that can tell them apart.

The map of the park, drawn from wait times alone

Here's the payoff. Plot every ride by its PC2 and PC3 loadings — its behavioral position — and color it by the park's category tag. Remember: the position comes purely from how each ride's queue moves over two seasons. The model was never told what any ride is or where it sits.

And yet the picture is the park. Kiddieland collapses into a tight cluster on the left. The evening thrill rides sink to the bottom-right. The daytime family flat-rides float to the top. The two blue dots stranded inside Kiddieland are the Family-tagged rides the queue data insists are really kiddie rides.

In [23]:

from matplotlib.patches import Ellipse

# the tight Kiddieland pile = rides loading hard-negative on PC2
kiddie_core = loadings.index[loadings["PC2"] < -0.13]
kx, ky = loadings.loc[kiddie_core, "PC2"], loadings.loc[kiddie_core, "PC3"]

fig, ax = plt.subplots(figsize=(11, 8.5))
cat_colors = {"Family": "tab:blue", "Kiddie": "tab:orange", "Thrill": "tab:red"}
for c, g in annotated.dropna(subset=["category"]).groupby("category"):
    ax.scatter(g["PC2"], g["PC3"], label=c, color=cat_colors.get(c, "grey"),
               s=80, alpha=0.75, edgecolor="white", zorder=3)

# one callout for the dense kiddie pile, instead of 13 overlapping labels
w, h = (kx.max() - kx.min()) + 0.07, (ky.max() - ky.min()) + 0.07
ax.add_patch(Ellipse((kx.mean(), ky.mean()), width=w, height=h, facecolor="tab:orange",
                     alpha=0.08, edgecolor="tab:orange", ls="--", lw=1.5, zorder=1))
ax.annotate("Kiddieland\n(tight cluster of kiddie rides)", (kx.mean(), ky.mean() - h / 2),
            ha="center", va="top", fontsize=9, color="tab:orange", fontweight="bold")

# label only the rides that carry the story
plain = {"giant_wheel": (6, 5), "grand_carousel": (8, -3), "impulse": (7, -3),
         "bumper_cars": (7, 4), "motor_boats": (7, 0), "paratrooper": (7, 0),
         "tea_cups": (7, 0), "tilt-a-whirl": (7, -2), "roto_jets": (-42, 6)}
for ride, (dx, dy) in plain.items():
    p = annotated.loc[ride]
    ax.annotate(ride, (p["PC2"], p["PC3"]), fontsize=8.5,
                xytext=(dx, dy), textcoords="offset points", zorder=4)

# the two Family-tagged rides that behave like kiddie rides — pull labels out with arrows
misfiled = {"s&g_carousel": (-70, -30), "kiddie_himalaya": (-78, 20)}
for ride, (dx, dy) in misfiled.items():
    p = annotated.loc[ride]
    ax.annotate(ride + "  (tagged Family)", (p["PC2"], p["PC3"]), fontsize=8.5,
                color="tab:blue", xytext=(dx, dy), textcoords="offset points",
                arrowprops=dict(arrowstyle="->", lw=0.7, color="tab:blue"), zorder=5)

ax.axhline(0, color="grey", lw=0.5)
ax.axvline(0, color="grey", lw=0.5)
ax.set_xlabel("PC2     Kiddieland corner  ←      →  main midway")
ax.set_ylabel("PC3     evening crowd  ←      →  daytime / school-group daypart")
ax.set_title("Knoebels rides in PC2 × PC3 space\n"
             "(position comes from queue patterns alone; color = the park's category tag)")
ax.legend(title="Category tag", loc="lower left")
plt.tight_layout()

What I took away

• PCA found real structure with zero domain input. Three axes — busyness, geography, time of day — fell out of nothing but wait times, and each one maps to something you could verify against the park independently (the average wait, the GPS map, the operating schedule).
• Standardizing matters twice. Once before PCA (so small rides aren't drowned out), and again before k-means on coordinates (so the clustering is reproducible).
• The model's labels beat the official ones. Where the category tag and the queue behavior disagreed (giant_wheel, s&g_carousel), the behavior was the truer description of how a ride is actually used.
• Be willing to kill a pretty hypothesis. The PC3 "throughput" story was tidy and wrong; the data pointed at time-of-day instead.

Skipping ahead in the textbook to chase a dataset I cared about turned out to be a great way to learn this.