kuh-no-bowls/knb_stats/pca_svd.ipynb

In [25]:

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"],
    )
)

# LZ_attractions_io carries one row per (timestamp, ride) and uses the
# Item _id that matches both ride_master and LZ_attractions_io_poi.
# (LZ_attractions_io_queuetimes is a same-data duplicate keyed by QueueLine
# id, which doesn't join to anything — don't use it.)
QUERY = text(
    """
    SELECT a.time_stamp,
           a._id            AS ride_id,
           r.ride           AS ride_name,
           r.capacity,
           r.duration,
           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
     ORDER BY a.time_stamp
    """
)

df = pd.read_sql(QUERY, engine)
df["queue_time_min"] = df["queue_time_sec"] / 60.0

# Scraper writes naive timestamps that represent America/New_York wall clock.
# Localize so downstream time-of-day / day-of-week features are unambiguous.
df["time_stamp"] = (
    pd.to_datetime(df["time_stamp"])
    .dt.tz_localize("America/New_York", ambiguous="NaT", nonexistent="NaT")
)
df = df.dropna(subset=["time_stamp"]).reset_index(drop=True)


def get_ride_time_matrix(df: pd.DataFrame, freq: str = "5min") -> pd.DataFrame:
    """Pivot the long df into a (timestamp x ride_name) matrix of queue times.

    Rows are snapped to `freq` bins so all rides share a common index.
    Values are mean queue time in minutes within each 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().fillna(0) # Impute missing values: forward fill first (assume line stayed the same), then fill remaining NaNs (like morning startup) with 0
    )


print(f"Loaded {len(df):,} rows across {df['ride_name'].nunique()} rides")
print(f"Time range: {df['time_stamp'].min()} → {df['time_stamp'].max()}")

Loaded 635,995 rows across 57 rides
Time range: 2024-05-15 12:23:34.954315-04:00 → 2026-05-26 18:00:20.956708-04:00

In [26]:

'''
if we were to run SVD on raw minutes, the algorithm will just chase the rides
with the biggest raw numbers.

Frist, need to standardize the data so every ride has a mean of 0 and a variance of 1.
scikit-learn makes this trivial
'''
from sklearn.preprocessing import StandardScaler

# Rides to exclude from the PCA:
#   powersurge — torn down years ago, ~90% of its rows are imputed zeros and
#   its loadings are pure artifact.
EXCLUDED_RIDES = ["powersurge"]

df_matrix = get_ride_time_matrix(df).drop(columns=EXCLUDED_RIDES, errors="ignore")

# Standardize the columns (rides)
scaler = StandardScaler()
scaled_matrix = scaler.fit_transform(df_matrix)

pd.DataFrame(scaled_matrix, index=df_matrix.index, columns=df_matrix.columns).head()

Out [26]:

ride_name	antique_cars	balloon_race	black_diamond	bumper_cars	cosmotron	downdraft	fandango	flyer	flying_tigers	flying_turns	...	spanish_bambini	stratosfear	super_round-up	tea_cups	tilt-a-whirl	tornado	tumbling_timbers	twister	umbrella_ride	whipper
ts
2024-05-15 12:20:00-04:00	-0.085045	1.45794	-0.063402	-0.214064	0.799599	1.160288	-2.167269	0.171589	2.954566	-0.970477	...	2.494992	0.589849	-5.171517	1.110184	-1.663817	1.064078	5.035897	-0.281509	1.255617	1.045173
2024-05-15 12:55:00-04:00	-0.085045	1.45794	-0.063402	-0.214064	0.799599	1.160288	-2.167269	0.171589	2.954566	-0.970477	...	2.494992	0.589849	-5.171517	1.110184	-1.663817	1.064078	5.035897	-0.281509	1.255617	1.045173
2024-05-15 13:00:00-04:00	-0.085045	1.45794	-0.063402	-0.214064	0.799599	1.160288	-2.167269	0.171589	2.954566	-0.970477	...	2.494992	0.589849	-5.171517	1.110184	-1.663817	1.064078	5.035897	-0.281509	1.255617	1.045173
2024-05-15 13:15:00-04:00	-0.085045	1.45794	-0.063402	-0.214064	0.799599	1.160288	-2.167269	0.171589	2.954566	-0.970477	...	-0.361261	0.589849	-5.171517	1.110184	-1.663817	1.064078	5.035897	-0.281509	1.255617	1.045173
2024-05-15 13:30:00-04:00	-0.085045	1.45794	-0.063402	-0.214064	0.799599	1.160288	-2.167269	0.171589	2.954566	-0.970477	...	-0.361261	0.589849	-5.171517	1.110184	-1.663817	1.064078	5.035897	-0.281509	1.255617	1.045173

5 rows × 56 columns

In [27]:

from sklearn.decomposition import PCA

# Let's extract the top 3 components ("Park Mood", "Capacity", etc.)
pca = PCA(n_components=3)
principal_components = pca.fit_transform(scaled_matrix)

# Look at how much variance (signal) each component explains
print(f"Variance explained by PC1 (Park Mood): {pca.explained_variance_ratio_[0]:.1%}")
print(f"Variance explained by PC2: {pca.explained_variance_ratio_[1]:.1%}")
print(f"Variance explained by PC3: {pca.explained_variance_ratio_[2]:.1%}")

Variance explained by PC1 (Park Mood): 40.4%
Variance explained by PC2: 13.0%
Variance explained by PC3: 6.7%

In [28]:

# Interpret each PC by looking at (1) which rides load onto it and (2) when its score is high/low.

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

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

# Top-magnitude rides per PC — these are the ones "shaping" the axis.
for pc in loadings.columns:
    print(f"\n=== {pc} — top 8 rides by |loading| ===")
    top = loadings[pc].reindex(loadings[pc].abs().sort_values(ascending=False).index).head(8)
    print(top.to_string())

# Daily-mean PC scores over time — reveals seasonality each PC captures.
scores.resample("D").mean().plot(subplots=True, figsize=(12, 6), title=["PC1", "PC2", "PC3"])

Out [28]:

=== PC1 — top 8 rides by |loading| ===
ride_name
stratosfear        0.179921
merry_mixer        0.174591
flyer              0.168501
s&g_carousel       0.167793
tornado            0.166206
kiddie_whip        0.165978
red_baron          0.165471
spanish_bambini    0.165410

=== PC2 — top 8 rides by |loading| ===
ride_name
bumper_cars        0.227223
giant_wheel        0.203075
kiddie_boats      -0.203019
s&g_carousel      -0.201848
spanish_bambini   -0.201826
pete's_fleet      -0.201538
helicopters       -0.197268
red_baron         -0.197263

=== PC3 — top 8 rides by |loading| ===
ride_name
paratrooper           0.343496
kiddie_bumper_cars   -0.328093
motor_boats          -0.268827
roto_jets             0.215735
stratosfear           0.196443
tornado               0.195098
impulse              -0.189184
tumbling_timbers     -0.188705

array([<Axes: title={'center': 'PC1'}, xlabel='ts'>,
       <Axes: title={'center': 'PC2'}, xlabel='ts'>,
       <Axes: title={'center': 'PC3'}, xlabel='ts'>], dtype=object)

In [29]:

# Join ride metadata (area + category) to the loadings so the labels write themselves.
# NOTE: attractions.io category names carry a trailing space ("Kiddie ", "Family ",
# "Thrill "), so trim() here — otherwise downstream == / .isin() / dict lookups on
# the bare names silently miss (e.g. the PC2×PC3 scatter colors all fall through).

metadata = pd.read_sql(
    text(
        """
        SELECT r.ride               AS ride_name,
               trim(c.name)         AS category,
               p.parent             AS area
          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")

annotated = loadings.join(metadata)

print("=== Mean loading by area ===")
print(annotated.groupby("area")[["PC1", "PC2", "PC3"]].mean().round(3))

print("\n=== Mean loading by category ===")
print(annotated.groupby("category")[["PC1", "PC2", "PC3"]].mean().round(3))

=== Mean loading by area ===
Empty DataFrame
Columns: [PC1, PC2, PC3]
Index: []

=== Mean loading by category ===
            PC1    PC2    PC3
category                     
Family    0.127  0.071 -0.016
Kiddie    0.135 -0.132 -0.080
Thrill    0.124  0.078  0.034

In [30]:

# Validate PC interpretations from the time side.
# PC1 (busyness) should peak weekend afternoons.
# PC2 (Kiddieland vs main midway) should fade in the evening as the kids leave.
# PC3 (daypart) should run positive midday — school/daycare group rates flood
# the daytime family rides — and go negative after dark, when the evening
# coaster/bumper-car/Ferris-wheel crowd takes over.
import matplotlib.pyplot as plt

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

scores.groupby(scores.index.hour).mean().plot(ax=axes[0])
axes[0].set_title("Mean PC score by hour of day")
axes[0].set_xlabel("Hour (local)")
axes[0].axhline(0, color="grey", linewidth=0.5)

scores.groupby(scores.index.dayofweek).mean().plot(ax=axes[1])
axes[1].set_title("Mean PC score by day of week (0=Mon, 6=Sun)")
axes[1].set_xlabel("Day of week")
axes[1].axhline(0, color="grey", linewidth=0.5)

plt.tight_layout()

In [31]:

# Low-rank reconstruction: how much of any one ride's queue is explained by the top 3 PCs?
# pca.inverse_transform takes us back to standardized space; scaler.inverse_transform back to minutes.
reconstructed = pd.DataFrame(
    scaler.inverse_transform(pca.inverse_transform(principal_components)),
    index=df_matrix.index,
    columns=df_matrix.columns,
)

# Per-ride R^2 of the rank-3 model: 1 - SS_residual / SS_total.
ss_resid = ((df_matrix - reconstructed) ** 2).sum()
ss_total = ((df_matrix - df_matrix.mean()) ** 2).sum()
ride_r2 = (1 - ss_resid / ss_total).rename("rank3_R2")

print("=== Best explained by top 3 PCs (move with the park) ===")
print(ride_r2.sort_values(ascending=False).head(10).round(3).to_string())
print("\n=== Worst explained (ride-specific dynamics) ===")
print(ride_r2.sort_values().head(10).round(3).to_string())

# Visualize one well-explained ride vs the reconstruction.
RIDE = ride_r2.sort_values(ascending=False).index[0]
fig, ax = plt.subplots(figsize=(14, 4))
df_matrix[RIDE].resample("D").mean().plot(ax=ax, label="actual", alpha=0.6)
reconstructed[RIDE].resample("D").mean().plot(ax=ax, label="top-3 PC reconstruction", alpha=0.9)
ax.set_title(f"{RIDE}: actual vs 3-PC reconstruction (daily-mean queue, min)")
ax.set_ylabel("Queue (min)")
ax.legend()

Out [31]:

=== Best explained by top 3 PCs (move with the park) ===
ride_name
s&g_carousel       0.945
spanish_bambini    0.936
kiddie_whip        0.913
red_baron          0.909
pete's_fleet       0.896
helicopters        0.882
stratosfear        0.879
pony_carts         0.862
kiddie_boats       0.858
goin'_buggy        0.842

=== Worst explained (ride-specific dynamics) ===
ride_name
paradrop             0.079
kiddie_firetrucks    0.197
ribbit               0.197
super_round-up       0.201
sklooosh             0.220
giant_flume          0.278
flying_tigers        0.279
scenic_skyway        0.286
ole_smokey           0.315
flying_turns         0.323

<matplotlib.legend.Legend at 0x7f62edf857f0>

In [32]:

# Blog-ready scatter: rides positioned in PC2 x PC3 space, colored by category.
# PC2 = Kiddieland corner vs main midway (spatial; see GPS-cluster validation)
# PC3 = evening crowd (−) vs daytime / school-group daypart (+); NOT throughput.
fig, ax = plt.subplots(figsize=(10, 8))

category_colors = {"Family": "tab:blue", "Kiddie": "tab:orange", "Thrill": "tab:red"}
for category, group in annotated.dropna(subset=["category"]).groupby("category"):
    ax.scatter(
        group["PC2"], group["PC3"],
        label=category,
        color=category_colors.get(category, "grey"),
        s=80, alpha=0.7, edgecolor="white",
    )

# Label only the rides farthest from the origin so the plot stays readable.
distance = (annotated["PC2"] ** 2 + annotated["PC3"] ** 2) ** 0.5
to_label = distance.sort_values(ascending=False).head(15).index
for ride in to_label:
    row = annotated.loc[ride]
    ax.annotate(ride, (row["PC2"], row["PC3"]), fontsize=8, alpha=0.85,
                xytext=(4, 2), textcoords="offset points")

ax.axhline(0, color="grey", linewidth=0.5)
ax.axvline(0, color="grey", linewidth=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 loading space")
ax.legend(title="Category")
plt.tight_layout()

In [33]:

# Derive park areas from GPS instead of a missing `parent` column.
# Each POI has a (lat, lon) point; kiddie rides are physically clustered in
# Kiddieland. KMeans on coords gives us an area label derived from geometry
# alone, which we can then cross-check against the PCA loadings.
#
# k=5 carves out a tight 19-ride Kiddieland. k=3 was too coarse — it swept
# main-midway rides near the park entrance (giant_wheel, grand_carousel,
# bumper_cars) into the kiddie blob. IMPORTANT: lat and lon must be standardized
# before KMeans — on raw degrees the clustering is seed-unstable (the contrast r
# wobbles 0.60-0.71 between random_states); standardizing makes it identical
# across seeds at r=0.78. The kiddie cluster is identified by lowest mean PC2
# loading, so this stays self-contained (no hand-typed ride list).
from sklearn.cluster import KMeans

poi = 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")

located = annotated.join(poi, how="inner")

# KMeans on the geometry only — no PCA information leaks in. Standardize the
# coordinates first so the clustering is reproducible (see note above).
N_AREAS = 5
coords = StandardScaler().fit_transform(located[["lat", "lon"]])
kmeans = KMeans(n_clusters=N_AREAS, n_init=10, random_state=0)
located["area_cluster"] = kmeans.fit_predict(coords)

# Identify the Kiddieland cluster as the one with the lowest mean PC2 loading,
# then export ride-name lists for the validation cell to build a contrast from
# geometry-derived groups (rather than the misfiling `category` tag).
kiddie_cluster = located.groupby("area_cluster")["PC2"].mean().idxmin()
kiddie_cluster_rides = located.index[located["area_cluster"] == kiddie_cluster]
main_cluster_rides   = located.index[located["area_cluster"] != kiddie_cluster]

fig, axes = plt.subplots(1, 2, figsize=(15, 7))

# Left: physical map colored by PC2 loading. If the kiddie hypothesis is right,
# one corner of the park should be uniformly blue.
sc = axes[0].scatter(
    located["lon"], located["lat"],
    c=located["PC2"], cmap="RdBu_r", s=140,
    edgecolor="black", linewidth=0.4,
)
axes[0].set_title("Rides colored by PC2 loading\n(blue = kiddie side, red = main park)")
axes[0].set_xlabel("Longitude")
axes[0].set_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.65,
                     xytext=(3, 3), textcoords="offset points")

# Right: same map, colored by KMeans cluster on coords alone.
for cluster, group in located.groupby("area_cluster"):
    label = f"Area {cluster}" + (" (Kiddieland)" if cluster == kiddie_cluster else "")
    axes[1].scatter(group["lon"], group["lat"], label=label, s=140, alpha=0.85)
axes[1].set_title(f"Spatial clusters (KMeans k={N_AREAS} on standardized lat/lon)")
axes[1].set_xlabel("Longitude")
axes[1].set_ylabel("Latitude")
axes[1].legend()

plt.tight_layout()

print(f"Kiddieland = cluster {kiddie_cluster} ({len(kiddie_cluster_rides)} rides):")
print("  " + ", ".join(sorted(kiddie_cluster_rides)))
print("\n=== Mean PC loading per spatial cluster ===")
print(located.groupby("area_cluster")[["PC1", "PC2", "PC3"]].mean().round(3))
print("\n=== Ride counts per cluster ===")
print(located.groupby("area_cluster").size())

Kiddieland = cluster 0 (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

=== Mean PC loading per spatial cluster ===
                PC1    PC2    PC3
area_cluster                     
0             0.138 -0.128 -0.075
1             0.116  0.057  0.094
2             0.129  0.078 -0.020
3             0.139  0.047  0.008
4             0.117  0.141 -0.063

=== Ride counts per cluster ===
area_cluster
0    19
1    12
2     6
3     8
4    11
dtype: int64

In [34]:

# Diagnostics: refine the PC labels, and check whether ffill+0 imputation is
# corrupting any rides into a "broken-ride" axis.

for pc in ["PC1", "PC2", "PC3"]:
    print(f"\n=== {pc} — full sorted loadings (most negative → most positive) ===")
    print(loadings[pc].sort_values().to_string())

silent_rate = (df_matrix == 0).mean().sort_values(ascending=False).rename("silent_rate")
print("\n=== Per-ride silent rate (fraction of rows = 0; high = closed often) ===")
print(silent_rate.head(20).round(3).to_string())

=== PC1 — full sorted loadings (most negative → most positive) ===
ride_name
paradrop              0.042544
giant_flume           0.044335
motor_boats           0.054094
sklooosh              0.055152
super_round-up        0.065455
kiddie_firetrucks     0.071166
flying_turns          0.072880
ribbit                0.079402
satellite             0.089203
scenic_skyway         0.090318
impulse               0.092897
flying_tigers         0.093864
tumbling_timbers      0.098942
ole_smokey            0.099403
kiddie_bumper_cars    0.099480
umbrella_ride         0.100910
hand_cars             0.104304
bumper_cars           0.106832
downdraft             0.107728
paratrooper           0.113176
giant_wheel           0.117456
pioneer_train         0.120164
antique_cars          0.123976
twister               0.135100
kozmo's_kurves        0.137001
cosmotron             0.139632
tea_cups              0.140882
tilt-a-whirl          0.141175
kiddie_wheel          0.142947
panther_cars          0.144494
phoenix               0.145333
grand_carousel        0.145648
looper                0.148817
balloon_race          0.149454
fandango              0.152827
kiddie_himalaya       0.154782
black_diamond         0.155503
galleon               0.155716
roto_jets             0.156365
kiddie_boats          0.156734
jet_skyfighter        0.157865
helicopters           0.159951
italian_trapeze       0.160205
pony_carts            0.160726
goin'_buggy           0.161893
whipper               0.162065
pete's_fleet          0.162164
haunted_mansion       0.162481
spanish_bambini       0.165410
red_baron             0.165471
kiddie_whip           0.165978
tornado               0.166206
s&g_carousel          0.167793
flyer                 0.168501
merry_mixer           0.174591
stratosfear           0.179921

=== PC2 — full sorted loadings (most negative → most positive) ===
ride_name
kiddie_boats         -0.203019
s&g_carousel         -0.201848
spanish_bambini      -0.201826
pete's_fleet         -0.201538
helicopters          -0.197268
red_baron            -0.197263
pony_carts           -0.194792
kiddie_whip          -0.194522
jet_skyfighter       -0.190656
goin'_buggy          -0.184707
kiddie_wheel         -0.180319
kiddie_himalaya      -0.149156
hand_cars            -0.131858
umbrella_ride        -0.130396
kiddie_firetrucks    -0.084228
kiddie_bumper_cars   -0.075835
paradrop             -0.071294
paratrooper          -0.071188
super_round-up       -0.042809
tumbling_timbers     -0.033110
panther_cars         -0.016679
roto_jets            -0.010444
tornado               0.000492
stratosfear           0.011899
fandango              0.015212
ole_smokey            0.029167
black_diamond         0.035973
phoenix               0.046732
ribbit                0.053584
balloon_race          0.059486
flyer                 0.069961
tea_cups              0.071883
downdraft             0.074681
haunted_mansion       0.074961
flying_tigers         0.079766
scenic_skyway         0.081828
italian_trapeze       0.083423
sklooosh              0.099405
looper                0.101226
merry_mixer           0.103398
antique_cars          0.115768
galleon               0.116481
flying_turns          0.119074
whipper               0.129719
kozmo's_kurves        0.138432
giant_flume           0.142256
twister               0.142868
pioneer_train         0.149265
tilt-a-whirl          0.159538
motor_boats           0.161983
satellite             0.174114
grand_carousel        0.175643
cosmotron             0.183703
impulse               0.185476
giant_wheel           0.203075
bumper_cars           0.227223

=== PC3 — full sorted loadings (most negative → most positive) ===
ride_name
kiddie_bumper_cars   -0.328093
motor_boats          -0.268827
impulse              -0.189184
tumbling_timbers     -0.188705
downdraft            -0.182160
giant_wheel          -0.175637
flying_turns         -0.162361
bumper_cars          -0.158399
super_round-up       -0.154787
ole_smokey           -0.150555
sklooosh             -0.144456
scenic_skyway        -0.118715
pioneer_train        -0.115972
kiddie_wheel         -0.108957
panther_cars         -0.108549
ribbit               -0.094803
flying_tigers        -0.093359
kiddie_firetrucks    -0.089730
hand_cars            -0.089480
kozmo's_kurves       -0.088087
umbrella_ride        -0.087754
spanish_bambini      -0.072432
helicopters          -0.071282
kiddie_whip          -0.061688
s&g_carousel         -0.053637
satellite            -0.049943
jet_skyfighter       -0.044744
red_baron            -0.040739
antique_cars         -0.040011
pete's_fleet         -0.035771
twister              -0.026674
kiddie_boats         -0.022203
grand_carousel       -0.018140
pony_carts           -0.011517
goin'_buggy          -0.006920
paradrop              0.019726
cosmotron             0.020171
fandango              0.041648
haunted_mansion       0.062773
galleon               0.065372
phoenix               0.066976
whipper               0.075824
balloon_race          0.081406
kiddie_himalaya       0.111578
black_diamond         0.119748
merry_mixer           0.121156
tea_cups              0.133413
tilt-a-whirl          0.144879
looper                0.145004
italian_trapeze       0.147921
giant_flume           0.150791
flyer                 0.151130
tornado               0.195098
stratosfear           0.196443
roto_jets             0.215735
paratrooper           0.343496

=== Per-ride silent rate (fraction of rows = 0; high = closed often) ===
ride_name
giant_flume          0.298
paratrooper          0.025
looper               0.018
kiddie_wheel         0.005
paradrop             0.004
fandango             0.004
pony_carts           0.004
kiddie_himalaya      0.002
tilt-a-whirl         0.002
super_round-up       0.002
flyer                0.000
giant_wheel          0.000
helicopters          0.000
jet_skyfighter       0.000
kiddie_firetrucks    0.000
flying_tigers        0.000
cosmotron            0.000
downdraft            0.000
black_diamond        0.000
bumper_cars          0.000

In [35]:

# Validate PC3 as a *time-of-day* (daypart) axis.
#
# The original "queue turnover / throughput" hypothesis is untestable here:
# ride_master.capacity and .duration are 100% NULL (0 of 60 rides), so there's
# no throughput to correlate against. And splitting rides by PC3 sign shows no
# queue-magnitude difference (neg-PC3 mean queue 9.8 min / busy-rate 0.37 vs
# pos-PC3 8.8 / 0.35) — PC3 is NOT "builds a line vs walk-on."
#
# What PC3 actually tracks is WHEN a ride is busy. The PC3 score runs strongly
# positive midday and negative after dark. Mechanism (park ops): discounted
# school/daycare group rates run from open to ~mid-afternoon and flood the
# daytime family flat-rides (tea_cups, tilt-a-whirl, roto_jets, italian_trapeze);
# the evening crowd (teens / date-night) shifts to coasters, bumper cars, motor
# boats and the Ferris wheel after dark (impulse, bumper_cars, giant_wheel).

afternoon_rides = loadings.index[loadings["PC3"] > 0.1]
evening_rides   = loadings.index[loadings["PC3"] < -0.1]
print(f"daytime rides (PC3 > +0.1): {len(afternoon_rides)}  -> {sorted(afternoon_rides)}")
print(f"evening rides (PC3 < -0.1): {len(evening_rides)}  -> {sorted(evening_rides)}")

hour = df_matrix.index.hour
# De-trend by the park-wide hourly mean so we see *shape* (when each group is
# busy relative to the whole park), not PC1's overall nightly ramp.
park_by_hour = df_matrix.mean(axis=1).groupby(hour).mean()

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

# Left: the PC3 score's own daypart profile — this IS the axis.
scores.groupby(scores.index.hour)["PC3"].mean().plot(ax=axes[0], marker="o")
axes[0].axhline(0, color="grey", linewidth=0.5)
axes[0].set_title("PC3 score by hour of day\n(+ = daytime pattern, − = evening pattern)")
axes[0].set_xlabel("Hour (local)")
axes[0].set_ylabel("Mean PC3 score")

# Right: each loading group's queue as a share of the park-wide level, by hour.
for name, grp, color in [("daytime rides (PC3 +)", afternoon_rides, "tab:orange"),
                         ("evening rides (PC3 −)", evening_rides, "tab:blue")]:
    share = df_matrix[grp].mean(axis=1).groupby(hour).mean() / park_by_hour
    share.plot(ax=axes[1], marker="o", label=name, color=color)
axes[1].axhline(1.0, color="grey", linewidth=0.5)
axes[1].set_title("Group queue ÷ park-wide queue, by hour\n(de-trended — isolates daypart shape from PC1)")
axes[1].set_xlabel("Hour (local)")
axes[1].set_ylabel("Share of park-wide mean queue")
axes[1].legend()

plt.tight_layout()

daytime rides (PC3 > +0.1): 13  -> ['black_diamond', 'flyer', 'giant_flume', 'italian_trapeze', 'kiddie_himalaya', 'looper', 'merry_mixer', 'paratrooper', 'roto_jets', 'stratosfear', 'tea_cups', 'tilt-a-whirl', 'tornado']
evening rides (PC3 < -0.1): 15  -> ['bumper_cars', 'downdraft', 'flying_turns', 'giant_wheel', 'impulse', 'kiddie_bumper_cars', 'kiddie_wheel', 'motor_boats', 'ole_smokey', 'panther_cars', 'pioneer_train', 'scenic_skyway', 'sklooosh', 'super_round-up', 'tumbling_timbers']

In [36]:

# Validate the PC labels against the most direct possible ground truth:
# the average queue time of the rides each PC is supposed to represent.
# PC1 should track park-wide mean queue.
# PC2 should track the Kiddieland-vs-main-midway contrast (NOT raw kiddie mean —
# on busy days both go up together, so raw kiddie mean is muddied by PC1).
#
# The kiddie/main split here comes from the GPS-cluster cell above
# (`kiddie_cluster_rides` / `main_cluster_rides`), NOT the `category` tag. The
# taxonomy misfiles ~3 rides across the divide (the big giant_wheel and
# grand_carousel are tagged Family but sit on the main midway; the small
# s&g_carousel sits in Kiddieland), which capped the category-based contrast at
# r=0.73. The geometry-derived grouping (k=5) lifts it to ~0.78 — and uses
# spatial information the PCA never saw, so it's an independent validation.
#
# Run the GPS-cluster cell first; this depends on its exported ride lists.
kiddie_cols = df_matrix.columns.intersection(kiddie_cluster_rides)
main_cols   = df_matrix.columns.intersection(main_cluster_rides)
print(f"Kiddieland rides in df_matrix: {len(kiddie_cols)}  |  main-park rides: {len(main_cols)}")

park_mean   = df_matrix.mean(axis=1)
kiddie_mean = df_matrix[kiddie_cols].mean(axis=1)
main_mean   = df_matrix[main_cols].mean(axis=1)
contrast    = main_mean - kiddie_mean

fig, axes = plt.subplots(1, 3, figsize=(18, 5))

axes[0].scatter(park_mean, scores["PC1"], s=4, alpha=0.25)
axes[0].set_xlabel("Park-wide mean queue (min)")
axes[0].set_ylabel("PC1 score")
axes[0].set_title(f"PC1 vs park-wide mean\nPearson r = {park_mean.corr(scores['PC1']):.3f}")

axes[1].scatter(kiddie_mean, scores["PC2"], s=4, alpha=0.25)
axes[1].set_xlabel("Kiddieland mean queue (min)")
axes[1].set_ylabel("PC2 score")
axes[1].set_title(f"PC2 vs raw Kiddieland mean\nPearson r = {kiddie_mean.corr(scores['PC2']):.3f}")

axes[2].scatter(contrast, scores["PC2"], s=4, alpha=0.25)
axes[2].set_xlabel("Main-midway mean − Kiddieland mean (min)")
axes[2].set_ylabel("PC2 score")
axes[2].set_title(f"PC2 vs main-vs-Kiddieland contrast\nPearson r = {contrast.corr(scores['PC2']):.3f}")

plt.tight_layout()

Kiddieland rides in df_matrix: 19  |  main-park rides: 37