wes/kuh-no-bowls
1
0 Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Add Streamlit app + knb_stats analysis notebooks

Browse Source 
Co-Authored-By: Claude Opus 4.7 <[email protected]>
This commit is contained in:
wesandClaude Opus 4.7 committed 2026-05-27 10:06:19 -04:00
1 parent dd6ff18b04
commit ade774f5aa
8 files changed
+1661 -1
Download Patch File Download Diff File Expand all files Collapse all files

No files matched your search

knb_stats/.gitignore
+8
View File
@@ -0,0 +1,8 @@
.env
.venv/
__pycache__/
*.pyc
.ipynb_checkpoints/
.ghost
post_assets/
preview.html
knb_stats/pca_svd.ipynb
+308
View File
@@ -0,0 +1,308 @@
{
"cells": [
{
"cell_type": "code",
"execution_count": 9,
"id": "cae65a79",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Loaded 635,281 rows across 57 rides\n",
"Time range: 2024-05-15 12:23:34.954315-04:00 → 2026-05-26 16:35:16.363800-04:00\n"
]
}
],
"source": [
"import os\n",
"from pathlib import Path\n",
"\n",
"import pandas as pd\n",
"from dotenv import load_dotenv\n",
"from sqlalchemy import create_engine, text\n",
"\n",
"load_dotenv(Path.cwd() / \".env\")\n",
"\n",
"engine = create_engine(\n",
" \"postgresql+psycopg2://{user}:{pw}@{host}:{port}/{db}\".format(\n",
" user=os.environ[\"PG_USER\"],\n",
" pw=os.environ[\"PG_PASSWORD\"],\n",
" host=os.environ[\"PG_HOST\"],\n",
" port=os.environ[\"PG_PORT\"],\n",
" db=os.environ[\"PG_DB\"],\n",
" )\n",
")\n",
"\n",
"# LZ_attractions_io carries one row per (timestamp, ride) and uses the\n",
"# Item _id that matches both ride_master and LZ_attractions_io_poi.\n",
"# (LZ_attractions_io_queuetimes is a same-data duplicate keyed by QueueLine\n",
"# id, which doesn't join to anything — don't use it.)\n",
"QUERY = text(\n",
" \"\"\"\n",
" SELECT a.time_stamp,\n",
" a._id AS ride_id,\n",
" r.ride AS ride_name,\n",
" r.capacity,\n",
" r.duration,\n",
" a.queue_time AS queue_time_sec\n",
" FROM knoebels.\"LZ_attractions_io\" a\n",
" JOIN knoebels.ride_master r\n",
" ON r.attractions_dot_io_id = a._id\n",
" WHERE a.queue_time IS NOT NULL\n",
" ORDER BY a.time_stamp\n",
" \"\"\"\n",
")\n",
"\n",
"df = pd.read_sql(QUERY, engine)\n",
"df[\"queue_time_min\"] = df[\"queue_time_sec\"] / 60.0\n",
"\n",
"# Scraper writes naive timestamps that represent America/New_York wall clock.\n",
"# Localize so downstream time-of-day / day-of-week features are unambiguous.\n",
"df[\"time_stamp\"] = (\n",
" pd.to_datetime(df[\"time_stamp\"])\n",
" .dt.tz_localize(\"America/New_York\", ambiguous=\"NaT\", nonexistent=\"NaT\")\n",
")\n",
"df = df.dropna(subset=[\"time_stamp\"]).reset_index(drop=True)\n",
"\n",
"\n",
"def get_ride_time_matrix(df: pd.DataFrame, freq: str = \"5min\") -> pd.DataFrame:\n",
" \"\"\"Pivot the long df into a (timestamp x ride_name) matrix of queue times.\n",
"\n",
" Rows are snapped to `freq` bins so all rides share a common index.\n",
" Values are mean queue time in minutes within each bin.\n",
" \"\"\"\n",
" snapped = df.assign(ts=df[\"time_stamp\"].dt.floor(freq))\n",
" return (\n",
" snapped.pivot_table(\n",
" index=\"ts\",\n",
" columns=\"ride_name\",\n",
" values=\"queue_time_min\",\n",
" aggfunc=\"mean\",\n",
" )\n",
" .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\n",
" )\n",
"\n",
"\n",
"print(f\"Loaded {len(df):,} rows across {df['ride_name'].nunique()} rides\")\n",
"print(f\"Time range: {df['time_stamp'].min()} → {df['time_stamp'].max()}\")\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "261c73fd",
"metadata": {},
"outputs": [],
"source": "'''\nif we were to run SVD on raw minutes, the algorithm will just chase the rides\nwith the biggest raw numbers.\n\nFrist, need to standardize the data so every ride has a mean of 0 and a variance of 1.\nscikit-learn makes this trivial\n'''\nfrom sklearn.preprocessing import StandardScaler\n\n# Rides to exclude from the PCA:\n# powersurge — torn down years ago, ~90% of its rows are imputed zeros and\n# its loadings are pure artifact.\nEXCLUDED_RIDES = [\"powersurge\"]\n\ndf_matrix = get_ride_time_matrix(df).drop(columns=EXCLUDED_RIDES, errors=\"ignore\")\n\n# Standardize the columns (rides)\nscaler = StandardScaler()\nscaled_matrix = scaler.fit_transform(df_matrix)\n\npd.DataFrame(scaled_matrix, index=df_matrix.index, columns=df_matrix.columns).head()\n"
},
{
"cell_type": "code",
"execution_count": 12,
"id": "ac1f62a3",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Variance explained by PC1 (Park Mood): 40.0%\n",
"Variance explained by PC2: 12.9%\n",
"Variance explained by PC3: 7.0%\n"
]
}
],
"source": [
"from sklearn.decomposition import PCA\n",
"\n",
"# Let's extract the top 3 components (\"Park Mood\", \"Capacity\", etc.)\n",
"pca = PCA(n_components=3)\n",
"principal_components = pca.fit_transform(scaled_matrix)\n",
"\n",
"# Look at how much variance (signal) each component explains\n",
"print(f\"Variance explained by PC1 (Park Mood): {pca.explained_variance_ratio_[0]:.1%}\")\n",
"print(f\"Variance explained by PC2: {pca.explained_variance_ratio_[1]:.1%}\")\n",
"print(f\"Variance explained by PC3: {pca.explained_variance_ratio_[2]:.1%}\")"
]
},
{
"cell_type": "code",
"execution_count": 14,
"id": "4663779d",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"\n",
"=== PC1 — top 8 rides by |loading| ===\n",
"ride_name\n",
"stratosfear 0.179701\n",
"merry_mixer 0.173788\n",
"flyer 0.168306\n",
"s&g_carousel 0.168079\n",
"tornado 0.166412\n",
"kiddie_whip 0.166132\n",
"spanish_bambini 0.165741\n",
"red_baron 0.165697\n",
"\n",
"=== PC2 — top 8 rides by |loading| ===\n",
"ride_name\n",
"bumper_cars 0.230985\n",
"giant_wheel 0.209181\n",
"kiddie_boats -0.196950\n",
"pete's_fleet -0.195781\n",
"spanish_bambini -0.195189\n",
"s&g_carousel -0.195019\n",
"red_baron -0.190284\n",
"helicopters -0.189507\n",
"\n",
"=== PC3 — top 8 rides by |loading| ===\n",
"ride_name\n",
"paratrooper 0.335514\n",
"kiddie_bumper_cars -0.327763\n",
"motor_boats -0.244509\n",
"powersurge 0.232955\n",
"tumbling_timbers -0.207141\n",
"roto_jets 0.199833\n",
"stratosfear 0.184504\n",
"tornado 0.180223\n"
]
},
{
"data": {
"text/plain": [
"array([<Axes: title={'center': 'PC1'}, xlabel='ts'>,\n",
" <Axes: title={'center': 'PC2'}, xlabel='ts'>,\n",
" <Axes: title={'center': 'PC3'}, xlabel='ts'>], dtype=object)"
]
},
"execution_count": 14,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAA9QAAAIzCAYAAADh4lEKAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjksIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvJkbTWQAAAAlwSFlzAAAPYQAAD2EBqD+naQAA1pVJREFUeJzs3QeYVOXVB/CzvbG7sLSlN6kqRRDECoqi2I2KGuwaG4kGNdHEHhRb+Exi1xjsXcSCiKCoKFVABem9lwW2Adu/57wz5847d+6d3vf/e559ts3O3J2Fufe8p7wpDQ0NDQQAAAAAAAAAAUkN7OYAAAAAAAAAgIAaAAAAAAAAIEjIUAMAAAAAAAAEAQE1AAAAAAAAQBAQUAMAAAAAAAAEAQE1AAAAAAAAQBAQUAMAAAAAAAAEAQE1AAAAAAAAQBAQUAMAAAAAAAAEAQE1AAAAAAAAQBAQUAMAAMS5SZMmUUpKivGWnZ1NPXr0oLFjx9LOnTvdbsuf33HHHdSrVy/Kzc2lvLw8GjhwII0fP572799v3G7+/Pl08803q+9lZGSo+wUAAIDApAd4ewAAAIiRhx56iLp06UKHDh2i2bNn03PPPUdTp06lpUuXquB5wYIFNGrUKKqoqKAxY8aoYJktXLiQHn30Ufruu+9o+vTp6mv8cy+//DL17duXunbtSqtWrcLfFQAAIEAIqAEAABLEGWecQYMGDVIfX3fdddS8eXOaOHEiTZkyRX3v/PPPp7S0NFq8eLHKUOsefvhheumll4zPb7rpJvrrX/9KOTk5KtONgBoAACBwKPkGAABIUCeffLJ6v379enrhhRdo69atKsA2B9OsdevWdM8997h9zsE0AAAABA8BNQAAQIJau3ates+Z6k8++UQFyBdeeGGsDwsAAKDRQEANAACQIEpLS2nPnj20ZcsWevfdd1VPNQfRZ511Fi1fvlwNKsvMzIz1YQIAADQa6KEGAABIECNGjHD7vFOnTvTmm29Su3btqKysjPLz82N2bAAAAI0RAmoAAIAE8cwzz6gsdHp6uuqB7tmzJ6WmOorNCgoKqLy8PNaHCAAA0KggoAYAAEgQgwcPNqZ8m/EgsiVLllB1dTXKvgEAAKIEPdQAAABJ4Oyzz6aDBw/Shx9+GOtDAQAAaDQQUAMAACSBG2+8kdq0aUO333675Z7Su3btovHjx8fk2AAAAJIVSr4BAACSQLNmzWjy5Mk0atQo6t+/P40ZM4YGDhyovrdo0SJ6++23aejQocbtN27cSK+//rr6eOHCheq9BNw87Ozyyy+Pye8BAACQSFIaGhoaYn0QAAAAYG/SpEl09dVX04IFC2x7qMX27dvpiSeeoM8//5w2bdqkhpb17t2bLrjgAho7dqwaXsZmzZpFw4cPt7yPk046SX0fAAAAvENADQAAAAAAABAE9FADAAAAAAAABAEBNQAAAAAAAEAQEFADAAAAAAAABAEBNQAAAAAAAEAQEFADAAAAAAAABAEBNQAAAAAAAEAQ0imO1dfX07Zt2yg/P59SUlJifTgAAAAAAACQ5BoaGqi8vJzatm1LqampiRtQczDdoUOHWB8GAAAAAAAANDKbN2+m9u3bJ25AzZlp+UUKCgpifTgAAAAAAACQ5MrKylRiV+LRhA2opcybg2kE1AAAAAAAABAt/rQdYygZAAAAAAAAQBAQUAPEicqqWlqwYS/V1zfE+lAAAAAAAMAPcV3yDdCYTPhiOb0xdxO9fMUgGtGndawPBwAAAAASUF1dHdXU1MT6MOJaRkYGpaWlheW+EFADxIkdpYfU++1ljvcAAAAAAIFs9bRjxw7av38/njQ/NG3alIqLi0PenhkBNUCcqHOWeqPkGwAAAAACJcF0q1atKDc3N+RAMZkXHg4cOEC7du1Sn7dp0yak+0NADRAn6pyt07XooQYAAACAQK4j6+qMYLp58+Z47nzIyclR7zmo5ucslPJvDCUDiKPVMoYMNQAAAAAEQnqmOTMN/pHnKtR+cwTUAHFW8o0MNQAAAAAEA2Xe0X+uEFADxFlAXVdfH+tDAQAAAAAAPyCgBogT9c6S7zrE0wAAAAAACQEBNUCckFlkyFADAAAAACQGBNQA8Vby7cxUAwAAAAAku6uuukr1M/NbZmYmHXbYYfTQQw9RbW2tMbj3xRdfpCFDhlCTJk3U/tGDBg2ip556Sm1/xZYtW0a/+93vqHPnzup++HvRgoAaIM5KvjGUDAAAAAAak9NPP522b99Oq1evpttvv50eeOABeuKJJ9T3Lr/8crrtttvo3HPPpW+++YaWLFlC9957L02ZMoWmT5+ubsOBddeuXenRRx+l4uLiqB479qEGiLMMNbbNAgAAAIBQcWb3YE1d1J/InIy0gCdoZ2VlGYHwTTfdRJMnT6ZPPvmEunXrRm+++SZ9/PHHKqAWnIk+55xzqKysTH1+9NFHqzd21113UTQhoAaIsx5qZKgBAAAAIFQcTPe578uoP5G/PTSScjNDCzNzcnKopKREBdM9e/Z0C6YFB+2FhYUUayj5BogTkpmWTLWV1TvLadbKXVE8KgAAAACA6GXVZ8yYQV9++SWdfPLJqgScA+p4hgw1QJyQYWTeAuqb31xEq3dV0A93nUztmuZE8egAAAAAIJFw6TVni2PxuIH67LPP1MCxmpoaqq+vp8suu0z1UfPX4x0CaoAEylCXVFar9/sqqxFQAwAAAIAtLokOtfQ6WoYPH07PPfecmvLdtm1bSk93HHePHj1oxYoVFM9Q8g0QZ1O+vQXUNXX1Pm8DAAAAAJBI8vLy1HZZHTt2NIJpxpnqVatWqYneVuXhpaWlFGsIqAESqOS7tg5bawEAAABA43DxxRfT6NGj6dJLL6VHHnmEFi5cSBs3blSl4CNGjFDbaLHq6mq1nRa/8cdbt25VH69Zsybix5gYNQAAjUB9ve8p37XOG0k2GwAAAAAgmcvW33rrLXrxxRfplVdeoYcfflhlsLt3705XXHEFjRzp6BHftm0bDRgwwPi5J598Ur2ddNJJNGvWrIgeIwJqgHgr+fYSLEuwjZJvAAAAAEgGkyZN8vr91NRUuvHGG9WbHd6XmkvAYwEl3wBxQoLkOmdZt9X35XUCATUAAAAAQOwhQw2QIBlqGUimboOhZAAAAG6++HU7Ldq0T23Zc+6AdtStZRM8QwAQcQioAeItQ20TLOu91QioAQAAXA5W19Ef315snCuP6tQMATUARAVKvgHihMTLtgE1MtQAAACWDtXUGcH074d0pA5FuXimACAqkKEGiBP1gWSoMeUbAADA8hw5/rwj1GRggMYoVoO5GvNzlZQZ6gPVtVRd6+o3BUgEEiTL1lh2e1Cr26KHGgAAwGMOSWqKY5sdgMYmIyNDvT9w4ECsDyVhyHMlz12wki5Dvbeymo579Gvq07aAPrzp2FgfDoDfJEi2iacxlAwAAMDHOTSNI2qARigtLY2aNm1Ku3btUp/n5uZicclLZpqDaX6u+Dnj5y4USRdQT1mylQ7W1NFPG/fF+lAAAiJVJ7YZagwlAwAAsISAGoCouLhYPQ0SVIN3HEzLcxaKpAuof91aGutDAAip5NtmG2oMJQMAAPCVoUa5NzRi3O7Qpk0batWqFdXU1MT6cOIal3mHmplO3oB6CwJqSPRts/zIUGPgBAAAgMd5MRUl3wAqUAxXsAiNbCgZT0levavC+BxT7iDRJnybh4/pMJQMAADA+3k0HQE1AERZUgXUa3e7gmmGQciQaNNJzR/rarTMNaZ8AwAAeFZxYSgZAERbUgXUS7e5l3sj6IBEoZdw66XdOmSoAQAAbM6jznNnKnqoASDKkiqg3lvp3nyPgBoShd42rZd/62rrkKEGAACwPI86F6ZR8g0A0ZZUAXXFoVq3z+22HwJIyAx1ve+ycAAAgMZIzpEYSgYA0ZZUAXVltXtAjXgaEoUeINtVVugLRHZBNwAAQGMk1V3ooQaAaEuqgLqiChlqSEx6mbddQF2jTf9GOwMAAIDFPtSY8g0AUZbUJd8IOiBR6P9WbTPUCKgBAAC8nkfTMJQMAKIsqQLqSlOGWu9LBYhn+r9Vu3+3esk3FosAAAA8z6O+MtRz15XQrrJDeOoAIGySu+Rby+gBxDM9hq6z+Xer/3vGUDIAAIDASr4Xb9pHl7w4l0b9ezY1IOkCAGGS1AE1sniQkCXffmSoMZQMAAAgsIB61srd6v2eiiqas64ETx8AhEVSBdQo+YZkCKhr/RhKZrdXNQAAQGM+j6Z66aHeUFJpfPzG3I1ROS4AiPz//XcXbKJvVzkWzGIhqQLqiqo6t8+RoYak2jarDhlqAAAAb+fRdC8Z6t+2lRkff7lsJ+0oRS81QKL73w/r6a8f/koPfrqMarRr5WhKpyRSUVXj9jl6qCFR6DE0B9Tc25ViWmXXM9dYLAIAAPA8R6baBNSHaupo7e4K9XHvNgV0VMem1ECo9gJIdKOP7kDPzVpL63ZX0jvzN9HlQzt73Gbp1lJ6fc5GlcUuO+SIF2fdOYxa5WeH5RiSJqDm7N2hGseqRHZGqvoYg5sgUZgDZP40zXRNoAfU+LcNAADg/7ZZK3eUq3Nr87xMmvqn4z0WrQEgMeVnZ9BtI7rTvVOW0VMzVtN5A9qpr0l8eNObi+ir33Z6/mBDHJR8f/fdd3T22WdT27Zt1YvSxx9/7PZ9zrDdd9991KZNG8rJyaERI0bQ6tWrKVIqq13l3oU5zicRfaaQIMwBslUGGiXfAAAA3s+jdkPJljnLvfu0LUAwDZBkLhnckbq2yKPC3AzarrVypKelUmZaKvHLwjn92tIb1w6h7/8yXL0V5WWG7fGDzlBXVlZSv3796JprrqELLrjA4/uPP/44/fvf/6ZXX32VunTpQvfeey+NHDmSfvvtN8rODk963WrCNz9p2Rlp6uM6bSoyQDwzB9BWATWGkgEAAFiTNj+7gPq37aVGQA0AySUjLZUmXT2Y2jTNVh/r/nJ6T7rLine truncated
"text/plain": [
"<Figure size 1200x600 with 3 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"# Interpret each PC by looking at (1) which rides load onto it and (2) when its score is high/low.\n",
"\n",
"loadings = pd.DataFrame(\n",
" pca.components_.T,\n",
" index=df_matrix.columns,\n",
" columns=[\"PC1\", \"PC2\", \"PC3\"],\n",
")\n",
"\n",
"scores = pd.DataFrame(\n",
" principal_components,\n",
" index=df_matrix.index,\n",
" columns=[\"PC1\", \"PC2\", \"PC3\"],\n",
")\n",
"\n",
"# Top-magnitude rides per PC — these are the ones \"shaping\" the axis.\n",
"for pc in loadings.columns:\n",
" print(f\"\\n=== {pc} — top 8 rides by |loading| ===\")\n",
" top = loadings[pc].reindex(loadings[pc].abs().sort_values(ascending=False).index).head(8)\n",
" print(top.to_string())\n",
"\n",
"# Daily-mean PC scores over time — reveals seasonality each PC captures.\n",
"scores.resample(\"D\").mean().plot(subplots=True, figsize=(12, 6), title=[\"PC1\", \"PC2\", \"PC3\"])\n"
]
},
{
"cell_type": "code",
"id": "95277639",
"source": "# Join ride metadata (area + category) to the loadings so the labels write themselves.\n# NOTE: attractions.io category names carry a trailing space (\"Kiddie \", \"Family \",\n# \"Thrill \"), so trim() here — otherwise downstream == / .isin() / dict lookups on\n# the bare names silently miss (e.g. the PC2×PC3 scatter colors all fall through).\n\nmetadata = pd.read_sql(\n text(\n \"\"\"\n SELECT r.ride AS ride_name,\n trim(c.name) AS category,\n p.parent AS area\n FROM knoebels.ride_master r\n JOIN knoebels.\"LZ_attractions_io_poi\" p\n ON p._id = r.attractions_dot_io_id\n LEFT JOIN knoebels.\"LZ_attractions_io_categories\" c\n ON c._id = p.category\n \"\"\"\n ),\n engine,\n).set_index(\"ride_name\")\n\nannotated = loadings.join(metadata)\n\nprint(\"=== Mean loading by area ===\")\nprint(annotated.groupby(\"area\")[[\"PC1\", \"PC2\", \"PC3\"]].mean().round(3))\n\nprint(\"\\n=== Mean loading by category ===\")\nprint(annotated.groupby(\"category\")[[\"PC1\", \"PC2\", \"PC3\"]].mean().round(3))\n",
"metadata": {},
"execution_count": null,
"outputs": []
},
{
"cell_type": "code",
"id": "94e2c1be",
"source": "# Validate PC interpretations from the time side.\n# PC1 (busyness) should peak weekend afternoons.\n# PC2 (Kiddieland vs main midway) should fade in the evening as the kids leave.\n# PC3 (daypart) should run positive midday — school/daycare group rates flood\n# the daytime family rides — and go negative after dark, when the evening\n# coaster/bumper-car/Ferris-wheel crowd takes over.\nimport matplotlib.pyplot as plt\n\nfig, axes = plt.subplots(1, 2, figsize=(14, 4))\n\nscores.groupby(scores.index.hour).mean().plot(ax=axes[0])\naxes[0].set_title(\"Mean PC score by hour of day\")\naxes[0].set_xlabel(\"Hour (local)\")\naxes[0].axhline(0, color=\"grey\", linewidth=0.5)\n\nscores.groupby(scores.index.dayofweek).mean().plot(ax=axes[1])\naxes[1].set_title(\"Mean PC score by day of week (0=Mon, 6=Sun)\")\naxes[1].set_xlabel(\"Day of week\")\naxes[1].axhline(0, color=\"grey\", linewidth=0.5)\n\nplt.tight_layout()\n",
"metadata": {},
"execution_count": null,
"outputs": []
},
{
"cell_type": "code",
"id": "8663040b",
"source": "# Low-rank reconstruction: how much of any one ride's queue is explained by the top 3 PCs?\n# pca.inverse_transform takes us back to standardized space; scaler.inverse_transform back to minutes.\nreconstructed = pd.DataFrame(\n scaler.inverse_transform(pca.inverse_transform(principal_components)),\n index=df_matrix.index,\n columns=df_matrix.columns,\n)\n\n# Per-ride R^2 of the rank-3 model: 1 - SS_residual / SS_total.\nss_resid = ((df_matrix - reconstructed) ** 2).sum()\nss_total = ((df_matrix - df_matrix.mean()) ** 2).sum()\nride_r2 = (1 - ss_resid / ss_total).rename(\"rank3_R2\")\n\nprint(\"=== Best explained by top 3 PCs (move with the park) ===\")\nprint(ride_r2.sort_values(ascending=False).head(10).round(3).to_string())\nprint(\"\\n=== Worst explained (ride-specific dynamics) ===\")\nprint(ride_r2.sort_values().head(10).round(3).to_string())\n\n# Visualize one well-explained ride vs the reconstruction.\nRIDE = ride_r2.sort_values(ascending=False).index[0]\nfig, ax = plt.subplots(figsize=(14, 4))\ndf_matrix[RIDE].resample(\"D\").mean().plot(ax=ax, label=\"actual\", alpha=0.6)\nreconstructed[RIDE].resample(\"D\").mean().plot(ax=ax, label=\"top-3 PC reconstruction\", alpha=0.9)\nax.set_title(f\"{RIDE}: actual vs 3-PC reconstruction (daily-mean queue, min)\")\nax.set_ylabel(\"Queue (min)\")\nax.legend()\n",
"metadata": {},
"execution_count": null,
"outputs": []
},
{
"cell_type": "code",
"id": "701ce1d6",
"source": "# Blog-ready scatter: rides positioned in PC2 x PC3 space, colored by category.\n# PC2 = Kiddieland corner vs main midway (spatial; see GPS-cluster validation)\n# PC3 = evening crowd (−) vs daytime / school-group daypart (+); NOT throughput.\nfig, ax = plt.subplots(figsize=(10, 8))\n\ncategory_colors = {\"Family\": \"tab:blue\", \"Kiddie\": \"tab:orange\", \"Thrill\": \"tab:red\"}\nfor category, group in annotated.dropna(subset=[\"category\"]).groupby(\"category\"):\n ax.scatter(\n group[\"PC2\"], group[\"PC3\"],\n label=category,\n color=category_colors.get(category, \"grey\"),\n s=80, alpha=0.7, edgecolor=\"white\",\n )\n\n# Label only the rides farthest from the origin so the plot stays readable.\ndistance = (annotated[\"PC2\"] ** 2 + annotated[\"PC3\"] ** 2) ** 0.5\nto_label = distance.sort_values(ascending=False).head(15).index\nfor ride in to_label:\n row = annotated.loc[ride]\n ax.annotate(ride, (row[\"PC2\"], row[\"PC3\"]), fontsize=8, alpha=0.85,\n xytext=(4, 2), textcoords=\"offset points\")\n\nax.axhline(0, color=\"grey\", linewidth=0.5)\nax.axvline(0, color=\"grey\", linewidth=0.5)\nax.set_xlabel(\"PC2 Kiddieland corner ← → main midway\")\nax.set_ylabel(\"PC3 evening crowd ← → daytime / school-group daypart\")\nax.set_title(\"Knoebels rides in PC2 × PC3 loading space\")\nax.legend(title=\"Category\")\nplt.tight_layout()\n",
"metadata": {},
"execution_count": null,
"outputs": []
},
{
"cell_type": "code",
"id": "6c2f885f",
"source": "# Derive park areas from GPS instead of a missing `parent` column.\n# Each POI has a (lat, lon) point; kiddie rides are physically clustered in\n# Kiddieland. KMeans on coords gives us an area label derived from geometry\n# alone, which we can then cross-check against the PCA loadings.\n#\n# k=5 carves out a tight 19-ride Kiddieland. k=3 was too coarse — it swept\n# main-midway rides near the park entrance (giant_wheel, grand_carousel,\n# bumper_cars) into the kiddie blob. IMPORTANT: lat and lon must be standardized\n# before KMeans — on raw degrees the clustering is seed-unstable (the contrast r\n# wobbles 0.60-0.71 between random_states); standardizing makes it identical\n# across seeds at r=0.78. The kiddie cluster is identified by lowest mean PC2\n# loading, so this stays self-contained (no hand-typed ride list).\nfrom sklearn.cluster import KMeans\n\npoi = pd.read_sql(\n text(\n \"\"\"\n SELECT r.ride AS ride_name,\n (p.location)[0]::float AS lat,\n (p.location)[1]::float AS lon\n FROM knoebels.ride_master r\n JOIN knoebels.\"LZ_attractions_io_poi\" p\n ON p._id = r.attractions_dot_io_id\n WHERE p.location IS NOT NULL\n \"\"\"\n ),\n engine,\n).set_index(\"ride_name\")\n\nlocated = annotated.join(poi, how=\"inner\")\n\n# KMeans on the geometry only — no PCA information leaks in. Standardize the\n# coordinates first so the clustering is reproducible (see note above).\nN_AREAS = 5\ncoords = StandardScaler().fit_transform(located[[\"lat\", \"lon\"]])\nkmeans = KMeans(n_clusters=N_AREAS, n_init=10, random_state=0)\nlocated[\"area_cluster\"] = kmeans.fit_predict(coords)\n\n# Identify the Kiddieland cluster as the one with the lowest mean PC2 loading,\n# then export ride-name lists for the validation cell to build a contrast from\n# geometry-derived groups (rather than the misfiling `category` tag).\nkiddie_cluster = located.groupby(\"area_cluster\")[\"PC2\"].mean().idxmin()\nkiddie_cluster_rides = located.index[located[\"area_cluster\"] == kiddie_cluster]\nmain_cluster_rides = located.index[located[\"area_cluster\"] != kiddie_cluster]\n\nfig, axes = plt.subplots(1, 2, figsize=(15, 7))\n\n# Left: physical map colored by PC2 loading. If the kiddie hypothesis is right,\n# one corner of the park should be uniformly blue.\nsc = axes[0].scatter(\n located[\"lon\"], located[\"lat\"],\n c=located[\"PC2\"], cmap=\"RdBu_r\", s=140,\n edgecolor=\"black\", linewidth=0.4,\n)\naxes[0].set_title(\"Rides colored by PC2 loading\\n(blue = kiddie side, red = main park)\")\naxes[0].set_xlabel(\"Longitude\")\naxes[0].set_ylabel(\"Latitude\")\nplt.colorbar(sc, ax=axes[0], label=\"PC2 loading\")\nfor ride, row in located.iterrows():\n axes[0].annotate(ride, (row[\"lon\"], row[\"lat\"]),\n fontsize=6, alpha=0.65,\n xytext=(3, 3), textcoords=\"offset points\")\n\n# Right: same map, colored by KMeans cluster on coords alone.\nfor cluster, group in located.groupby(\"area_cluster\"):\n label = f\"Area {cluster}\" + (\" (Kiddieland)\" if cluster == kiddie_cluster else \"\")\n axes[1].scatter(group[\"lon\"], group[\"lat\"], label=label, s=140, alpha=0.85)\naxes[1].set_title(f\"Spatial clusters (KMeans k={N_AREAS} on standardized lat/lon)\")\naxes[1].set_xlabel(\"Longitude\")\naxes[1].set_ylabel(\"Latitude\")\naxes[1].legend()\n\nplt.tight_layout()\n\nprint(f\"Kiddieland = cluster {kiddie_cluster} ({len(kiddie_cluster_rides)} rides):\")\nprint(\" \" + \", \".join(sorted(kiddie_cluster_rides)))\nprint(\"\\n=== Mean PC loading per spatial cluster ===\")\nprint(located.groupby(\"area_cluster\")[[\"PC1\", \"PC2\", \"PC3\"]].mean().round(3))\nprint(\"\\n=== Ride counts per cluster ===\")\nprint(located.groupby(\"area_cluster\").size())\n",
"metadata": {},
"execution_count": null,
"outputs": []
},
{
"cell_type": "code",
"id": "6c723cc6",
"source": "# Diagnostics: refine the PC labels, and check whether ffill+0 imputation is\n# corrupting any rides into a \"broken-ride\" axis.\n\nfor pc in [\"PC1\", \"PC2\", \"PC3\"]:\n print(f\"\\n=== {pc} — full sorted loadings (most negative → most positive) ===\")\n print(loadings[pc].sort_values().to_string())\n\nsilent_rate = (df_matrix == 0).mean().sort_values(ascending=False).rename(\"silent_rate\")\nprint(\"\\n=== Per-ride silent rate (fraction of rows = 0; high = closed often) ===\")\nprint(silent_rate.head(20).round(3).to_string())\n",
"metadata": {},
"execution_count": null,
"outputs": []
},
{
"cell_type": "code",
"id": "4c567eec",
"source": "# Validate PC3 as a *time-of-day* (daypart) axis.\n#\n# The original \"queue turnover / throughput\" hypothesis is untestable here:\n# ride_master.capacity and .duration are 100% NULL (0 of 60 rides), so there's\n# no throughput to correlate against. And splitting rides by PC3 sign shows no\n# queue-magnitude difference (neg-PC3 mean queue 9.8 min / busy-rate 0.37 vs\n# pos-PC3 8.8 / 0.35) — PC3 is NOT \"builds a line vs walk-on.\"\n#\n# What PC3 actually tracks is WHEN a ride is busy. The PC3 score runs strongly\n# positive midday and negative after dark. Mechanism (park ops): discounted\n# school/daycare group rates run from open to ~mid-afternoon and flood the\n# daytime family flat-rides (tea_cups, tilt-a-whirl, roto_jets, italian_trapeze);\n# the evening crowd (teens / date-night) shifts to coasters, bumper cars, motor\n# boats and the Ferris wheel after dark (impulse, bumper_cars, giant_wheel).\n\nafternoon_rides = loadings.index[loadings[\"PC3\"] > 0.1]\nevening_rides = loadings.index[loadings[\"PC3\"] < -0.1]\nprint(f\"daytime rides (PC3 > +0.1): {len(afternoon_rides)} -> {sorted(afternoon_rides)}\")\nprint(f\"evening rides (PC3 < -0.1): {len(evening_rides)} -> {sorted(evening_rides)}\")\n\nhour = df_matrix.index.hour\n# De-trend by the park-wide hourly mean so we see *shape* (when each group is\n# busy relative to the whole park), not PC1's overall nightly ramp.\npark_by_hour = df_matrix.mean(axis=1).groupby(hour).mean()\n\nfig, axes = plt.subplots(1, 2, figsize=(15, 5))\n\n# Left: the PC3 score's own daypart profile — this IS the axis.\nscores.groupby(scores.index.hour)[\"PC3\"].mean().plot(ax=axes[0], marker=\"o\")\naxes[0].axhline(0, color=\"grey\", linewidth=0.5)\naxes[0].set_title(\"PC3 score by hour of day\\n(+ = daytime pattern, − = evening pattern)\")\naxes[0].set_xlabel(\"Hour (local)\")\naxes[0].set_ylabel(\"Mean PC3 score\")\n\n# Right: each loading group's queue as a share of the park-wide level, by hour.\nfor name, grp, color in [(\"daytime rides (PC3 +)\", afternoon_rides, \"tab:orange\"),\n (\"evening rides (PC3 −)\", evening_rides, \"tab:blue\")]:\n share = df_matrix[grp].mean(axis=1).groupby(hour).mean() / park_by_hour\n share.plot(ax=axes[1], marker=\"o\", label=name, color=color)\naxes[1].axhline(1.0, color=\"grey\", linewidth=0.5)\naxes[1].set_title(\"Group queue ÷ park-wide queue, by hour\\n(de-trended — isolates daypart shape from PC1)\")\naxes[1].set_xlabel(\"Hour (local)\")\naxes[1].set_ylabel(\"Share of park-wide mean queue\")\naxes[1].legend()\n\nplt.tight_layout()\n",
"metadata": {},
"execution_count": null,
"outputs": []
},
{
"cell_type": "code",
"id": "6d7970ab",
"source": "# Validate the PC labels against the most direct possible ground truth:\n# the average queue time of the rides each PC is supposed to represent.\n# PC1 should track park-wide mean queue.\n# PC2 should track the Kiddieland-vs-main-midway contrast (NOT raw kiddie mean —\n# on busy days both go up together, so raw kiddie mean is muddied by PC1).\n#\n# The kiddie/main split here comes from the GPS-cluster cell above\n# (`kiddie_cluster_rides` / `main_cluster_rides`), NOT the `category` tag. The\n# taxonomy misfiles ~3 rides across the divide (the big giant_wheel and\n# grand_carousel are tagged Family but sit on the main midway; the small\n# s&g_carousel sits in Kiddieland), which capped the category-based contrast at\n# r=0.73. The geometry-derived grouping (k=5) lifts it to ~0.78 — and uses\n# spatial information the PCA never saw, so it's an independent validation.\n#\n# Run the GPS-cluster cell first; this depends on its exported ride lists.\nkiddie_cols = df_matrix.columns.intersection(kiddie_cluster_rides)\nmain_cols = df_matrix.columns.intersection(main_cluster_rides)\nprint(f\"Kiddieland rides in df_matrix: {len(kiddie_cols)} | main-park rides: {len(main_cols)}\")\n\npark_mean = df_matrix.mean(axis=1)\nkiddie_mean = df_matrix[kiddie_cols].mean(axis=1)\nmain_mean = df_matrix[main_cols].mean(axis=1)\ncontrast = main_mean - kiddie_mean\n\nfig, axes = plt.subplots(1, 3, figsize=(18, 5))\n\naxes[0].scatter(park_mean, scores[\"PC1\"], s=4, alpha=0.25)\naxes[0].set_xlabel(\"Park-wide mean queue (min)\")\naxes[0].set_ylabel(\"PC1 score\")\naxes[0].set_title(f\"PC1 vs park-wide mean\\nPearson r = {park_mean.corr(scores['PC1']):.3f}\")\n\naxes[1].scatter(kiddie_mean, scores[\"PC2\"], s=4, alpha=0.25)\naxes[1].set_xlabel(\"Kiddieland mean queue (min)\")\naxes[1].set_ylabel(\"PC2 score\")\naxes[1].set_title(f\"PC2 vs raw Kiddieland mean\\nPearson r = {kiddie_mean.corr(scores['PC2']):.3f}\")\n\naxes[2].scatter(contrast, scores[\"PC2\"], s=4, alpha=0.25)\naxes[2].set_xlabel(\"Main-midway mean − Kiddieland mean (min)\")\naxes[2].set_ylabel(\"PC2 score\")\naxes[2].set_title(f\"PC2 vs main-vs-Kiddieland contrast\\nPearson r = {contrast.corr(scores['PC2']):.3f}\")\n\nplt.tight_layout()\n",
"metadata": {},
"execution_count": null,
"outputs": []
}
],
"metadata": {
"kernelspec": {
"display_name": ".venv",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.14.4"
}
},
"nbformat": 4,
"nbformat_minor": 5
}
knb_stats/writeup.ipynb
+836
View File
@@ -0,0 +1,836 @@
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Reading a Theme Park's Mind with PCA\n",
"\n",
"### What 600,000 queue-time readings can tell you about how people move through Knoebels\n",
"\n",
"I'm working my way through *An Introduction to Statistical Learning* and Gilbert\n",
"Strang's linear algebra lectures, and I got impatient. Principal Component\n",
"Analysis sits a few chapters ahead of where I am, but I had a dataset I couldn't\n",
"stop thinking about, so I skipped ahead and used it as an excuse to learn PCA by\n",
"*doing*. This notebook is me figuring it out live.\n",
"\n",
"**The data.** Every few minutes, a scraper records the posted wait time for each\n",
"ride at [Knoebels](https://knoebels.com) (a wonderful, weird old amusement park\n",
"in Pennsylvania). Two seasons of that adds up to ~600k readings across ~57 rides.\n",
"\n",
"**The question.** If I *don't* tell the algorithm anything about the rides \u2014 not\n",
"which are roller coasters, not which are for toddlers, not where they sit in the\n",
"park \u2014 can it still discover the hidden \"axes\" that explain how the whole park\n",
"behaves? That's exactly what PCA is for: it finds the handful of directions that\n",
"explain the most variation in the data, and leaves us to interpret them.\n",
"\n",
"Spoiler: it finds three, and they turn out to be **how busy the park is**,\n",
"**which half of the park you're in**, and **what time of day it is** \u2014 none of\n",
"which it was ever told."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## The data: one big table of waits\n",
"\n",
"The natural shape for this is a table where every **row is a moment in time** (a\n",
"5-minute snapshot of the park) and every **column is a ride**. Each cell is the\n",
"wait, in minutes, for that ride at that moment.\n",
"\n",
"A couple of small judgment calls in building it:\n",
"\n",
"- **Forward-fill, then zero.** If a ride didn't report at 2:05 we carry its 2:00\n",
" value forward (the line probably didn't teleport). Anything still missing \u2014\n",
" like early morning before the ride opens \u2014 becomes 0.\n",
"- **Drop `powersurge`.** It hasn't been running for most of the time since I began collecting data; ~90% of its column is\n",
" imputed zeros, so it only adds noise."
]
},
{
"cell_type": "code",
"execution_count": 16,
"id": "80e904a5",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"635,995 readings -> matrix of 12,073 time-snapshots x 56 rides\n"
]
},
{
"data": {
"text/html": [
"<div>\n",
"<style scoped>\n",
" .dataframe tbody tr th:only-of-type {\n",
" vertical-align: middle;\n",
" }\n",
"\n",
" .dataframe tbody tr th {\n",
" vertical-align: top;\n",
" }\n",
"\n",
" .dataframe thead th {\n",
" text-align: right;\n",
" }\n",
"</style>\n",
"<table border=\"1\" class=\"dataframe\">\n",
" <thead>\n",
" <tr style=\"text-align: right;\">\n",
" <th>ride_name</th>\n",
" <th>antique_cars</th>\n",
" <th>balloon_race</th>\n",
" <th>black_diamond</th>\n",
" <th>bumper_cars</th>\n",
" <th>cosmotron</th>\n",
" </tr>\n",
" <tr>\n",
" <th>ts</th>\n",
" <th></th>\n",
" <th></th>\n",
" <th></th>\n",
" <th></th>\n",
" <th></th>\n",
" </tr>\n",
" </thead>\n",
" <tbody>\n",
" <tr>\n",
" <th>2024-05-15 12:20:00</th>\n",
" <td>10.0</td>\n",
" <td>10.0</td>\n",
" <td>10.0</td>\n",
" <td>10.0</td>\n",
" <td>10.0</td>\n",
" </tr>\n",
" <tr>\n",
" <th>2024-05-15 12:55:00</th>\n",
" <td>10.0</td>\n",
" <td>10.0</td>\n",
" <td>10.0</td>\n",
" <td>10.0</td>\n",
" <td>10.0</td>\n",
" </tr>\n",
" <tr>\n",
" <th>2024-05-15 13:00:00</th>\n",
" <td>10.0</td>\n",
" <td>10.0</td>\n",
" <td>10.0</td>\n",
" <td>10.0</td>\n",
" <td>10.0</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
"</div>"
],
"text/plain": [
"ride_name antique_cars balloon_race black_diamond bumper_cars \\\n",
"ts \n",
"2024-05-15 12:20:00 10.0 10.0 10.0 10.0 \n",
"2024-05-15 12:55:00 10.0 10.0 10.0 10.0 \n",
"2024-05-15 13:00:00 10.0 10.0 10.0 10.0 \n",
"\n",
"ride_name cosmotron \n",
"ts \n",
"2024-05-15 12:20:00 10.0 \n",
"2024-05-15 12:55:00 10.0 \n",
"2024-05-15 13:00:00 10.0 "
]
},
"execution_count": 16,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"import os\n",
"from pathlib import Path\n",
"\n",
"import pandas as pd\n",
"from dotenv import load_dotenv\n",
"from sqlalchemy import create_engine, text\n",
"\n",
"load_dotenv(Path.cwd() / \".env\")\n",
"engine = create_engine(\n",
" \"postgresql+psycopg2://{user}:{pw}@{host}:{port}/{db}\".format(\n",
" user=os.environ[\"PG_USER\"], pw=os.environ[\"PG_PASSWORD\"],\n",
" host=os.environ[\"PG_HOST\"], port=os.environ[\"PG_PORT\"], db=os.environ[\"PG_DB\"],\n",
" )\n",
")\n",
"\n",
"df = pd.read_sql(\n",
" text(\n",
" \"\"\"\n",
" SELECT a.time_stamp,\n",
" r.ride AS ride_name,\n",
" a.queue_time AS queue_time_sec\n",
" FROM knoebels.\"LZ_attractions_io\" a\n",
" JOIN knoebels.ride_master r ON r.attractions_dot_io_id = a._id\n",
" WHERE a.queue_time IS NOT NULL\n",
" \"\"\"\n",
" ),\n",
" engine,\n",
")\n",
"df[\"queue_time_min\"] = df[\"queue_time_sec\"] / 60.0\n",
"df[\"time_stamp\"] = pd.to_datetime(df[\"time_stamp\"])\n",
"\n",
"\n",
"def ride_time_matrix(df, freq=\"5min\"):\n",
" \"\"\"Long table -> (time x ride) matrix of mean wait in minutes per time bin.\"\"\"\n",
" snapped = df.assign(ts=df[\"time_stamp\"].dt.floor(freq))\n",
" return (\n",
" snapped.pivot_table(index=\"ts\", columns=\"ride_name\",\n",
" values=\"queue_time_min\", aggfunc=\"mean\")\n",
" .sort_index()\n",
" .ffill() # carry the last known wait forward\n",
" .fillna(0) # anything still missing (pre-open) = no line\n",
" )\n",
"\n",
"\n",
"df_matrix = ride_time_matrix(df).drop(columns=[\"powersurge\"], errors=\"ignore\")\n",
"print(f\"{len(df):,} readings -> matrix of {df_matrix.shape[0]:,} time-snapshots x {df_matrix.shape[1]} rides\")\n",
"df_matrix.iloc[:3, :5]"
]
},
{
"cell_type": "markdown",
"id": "89de2f88",
"metadata": {},
"source": [
"## Step 1: put every ride on the same footing\n",
"\n",
"Here's the first thing that bit me. PCA hunts for directions of **maximum\n",
"variance** \u2014 but \"variance\" is measured in whatever units the data is in. A\n",
"headline roller coaster might swing between 5 and 45 minutes; a kiddie ride\n",
"swings between 0 and 6. If I feed in raw minutes, PCA will conclude the coaster\n",
"is \"more important\" purely because its *numbers* are bigger, and the small rides\n",
"will be invisible.\n",
"\n",
"The fix is to **standardize** every column: subtract its mean and divide by its\n",
"standard deviation, so each ride has mean 0 and variance 1. Now a \"big day\" for\n",
"the kiddie carousel counts just as much as a big day for the coaster. We're\n",
"asking about *patterns of co-movement*, not raw magnitudes."
]
},
{
"cell_type": "code",
"execution_count": 17,
"id": "e3c1738c",
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"<div>\n",
"<style scoped>\n",
" .dataframe tbody tr th:only-of-type {\n",
" vertical-align: middle;\n",
" }\n",
"\n",
" .dataframe tbody tr th {\n",
" vertical-align: top;\n",
" }\n",
"\n",
" .dataframe thead th {\n",
" text-align: right;\n",
" }\n",
"</style>\n",
"<table border=\"1\" class=\"dataframe\">\n",
" <thead>\n",
" <tr style=\"text-align: right;\">\n",
" <th>ride_name</th>\n",
" <th>antique_cars</th>\n",
" <th>balloon_race</th>\n",
" <th>black_diamond</th>\n",
" <th>bumper_cars</th>\n",
" </tr>\n",
" </thead>\n",
" <tbody>\n",
" <tr>\n",
" <th>mean</th>\n",
" <td>0.0</td>\n",
" <td>-0.0</td>\n",
" <td>-0.0</td>\n",
" <td>-0.0</td>\n",
" </tr>\n",
" <tr>\n",
" <th>std</th>\n",
" <td>1.0</td>\n",
" <td>1.0</td>\n",
" <td>1.0</td>\n",
" <td>1.0</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
"</div>"
],
"text/plain": [
"ride_name antique_cars balloon_race black_diamond bumper_cars\n",
"mean 0.0 -0.0 -0.0 -0.0\n",
"std 1.0 1.0 1.0 1.0"
]
},
"execution_count": 17,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"from sklearn.preprocessing import StandardScaler\n",
"\n",
"scaler = StandardScaler()\n",
"scaled = scaler.fit_transform(df_matrix) # each column now mean 0, variance 1\n",
"pd.DataFrame(scaled, columns=df_matrix.columns).describe().loc[[\"mean\", \"std\"]].iloc[:, :4].round(2)"
]
},
{
"cell_type": "markdown",
"id": "dc6ae7d7",
"metadata": {},
"source": [
"## Step 2: ask for the axes of variation\n",
"\n",
"This is the part I'd been intimidated by, and it turned out to be a clean idea.\n",
"\n",
"We have ~57 ride-columns, but they're hugely redundant \u2014 when the park is\n",
"mobbed, *everything* gets busy together. PCA asks: **can I describe most of that\n",
"57-dimensional wobble with just a few new axes?** Each axis (a *principal\n",
"component*) is a weighted blend of the original rides, chosen so the first one\n",
"captures as much variance as possible, the second captures as much of what's\n",
"left as possible while being perpendicular to the first, and so on.\n",
"\n",
"The linear algebra underneath is the **singular value decomposition**. Our\n",
"standardized matrix $X$ factors as\n",
"\n",
"$$X = U\\,\\Sigma\\,V^{\\top}.$$\n",
"\n",
"The columns of $V$ are the new axes \u2014 the **loadings**, telling us how much each\n",
"ride contributes to each component. $U\\Sigma$ gives the **scores** \u2014 where each\n",
"time-snapshot lands along those axes. And the singular values in $\\Sigma$, squared,\n",
"tell us how much variance each axis explains. `sklearn`'s `PCA` does exactly this\n",
"and hands back the pieces with friendly names.\n",
"\n",
"I asked for the top 3 components. Look at how concentrated the signal is:"
]
},
{
"cell_type": "code",
"execution_count": 18,
"id": "404f4039",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"PC1: 40.4% of all variance\n",
"PC2: 13.0% of all variance\n",
"PC3: 6.7% of all variance\n",
" 60.1% captured by just 3 axes (out of 56)\n"
]
}
],
"source": [
"from sklearn.decomposition import PCA\n",
"\n",
"pca = PCA(n_components=3)\n",
"principal_components = pca.fit_transform(scaled)\n",
"\n",
"scores = pd.DataFrame(principal_components, index=df_matrix.index,\n",
" columns=[\"PC1\", \"PC2\", \"PC3\"])\n",
"loadings = pd.DataFrame(pca.components_.T, index=df_matrix.columns,\n",
" columns=[\"PC1\", \"PC2\", \"PC3\"])\n",
"\n",
"for i, pc in enumerate([\"PC1\", \"PC2\", \"PC3\"]):\n",
" print(f\"{pc}: {pca.explained_variance_ratio_[i]:5.1%} of all variance\")\n",
"print(f\"{'':4} {pca.explained_variance_ratio_[:3].sum():5.1%} captured by just 3 axes (out of 56)\")"
]
},
{
"cell_type": "markdown",
"id": "c8f88158",
"metadata": {},
"source": [
"Three axes out of fifty-six explain ~60% of everything. The rest is\n",
"ride-specific noise. Now the fun part \u2014 **what are these three axes?** PCA won't\n",
"tell us; it just hands us numbers. We have to read them.\n",
"\n",
"My approach for each component: look at **which rides load heavily on it** (the\n",
"`loadings`), form a hypothesis about what real-world thing that represents, then\n",
"**validate** the hypothesis against something the PCA never saw."
]
},
{
"cell_type": "markdown",
"id": "ef23e2c7",
"metadata": {},
"source": [
"## PC1 \u2014 the whole park breathes together\n",
"\n",
"The first component is the easiest. Every single ride loads on it with the\n",
"**same sign** \u2014 they all move together. That's the signature of one global force:\n",
"overall busyness. Empty Tuesday morning, the whole park is near zero; packed\n",
"Saturday afternoon, everything spikes at once.\n",
"\n",
"**The test:** if PC1 is really \"busyness,\" then a snapshot's PC1 score should\n",
"track the plain average wait across all rides at that moment. To make that easy\n",
"to eyeball, I put both quantities on the same standardized scale \u2014 so **if PC1\n",
"were exactly the average wait, every dot would land on the dashed 45\u00b0 line.** The\n",
"red line is the actual relationship (its slope is the correlation). For PC1 the\n",
"red line sits almost on top of the diagonal:"
]
},
{
"cell_type": "code",
"execution_count": 19,
"id": "72b3033b",
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAnQAAAKACAYAAAAVauXAAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjksIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvJkbTWQAAAAlwSFlzAAAPYQAAD2EBqD+naQAA5QRJREFUeJzsnQd4U1Ubx/9tk3Q33RMoUPaeIoIgDkAQFAUUByA4AUXAgYqTT1FURBEcyHIxBNmKKKBMWbJnKaXQvZvuJmm+5z31pkmajnSmyft7nvT2nNzcnNysf97poNPpdGAYhmEYhmEaLY4NvQCGYRiGYRimZrCgYxiGYRiGaeSwoGMYhmEYhmnksKBjGIZhGIZp5LCgYxiGYRiGaeSwoGMYhmEYhmnksKBjGIZhGIZp5LCgYxiGYRiGaeSwoGMYhmEYhmnksKBjmFrir7/+goODg9haG9a8tqqycuVK8RiOHTtWq8dt3rw57rnnHtgS9JgmTpwIa+a2224TF3t47VYFe3mcTN3Bgo6xqi9r6eLi4oI2bdpg2rRpSEpKKrM/zb344oto164d3Nzc4O7ujp49e+J///sfMjMz9fsdOXIEU6ZMEdfJ5XJxbKaEn376CQsXLqzy6Xj//fexadOmRnn6zp8/j7fffhvXrl1r6KUwds6SJUvE511dvEdNodf7888/j7Zt24rPyfbt2+Orr76q9vEY60bW0AtgGEPeffddtGjRAgUFBdi/fz++/PJL/Prrrzh79qz4QCKOHj2KYcOGIScnB48++qgQawRZbj744APs3bsXO3fuFHN022+//RZdunRBy5YtcfnyZbs84QMGDEB+fj4UCoXRlwWd1xdeeKHKgm706NG477770BgF3TvvvCMsQmS9YpiGFHT+/v5lLKi18R41hX7E/PPPP3jssccQGBiI9evX49lnnxU/mK3dgstYDgs6xqq4++670atXL/H/E088AT8/PyxYsACbN2/GuHHjhPVt1KhRcHJywokTJ4SFzpD33nsPS5cu1Y/pw+uVV16Bq6ursPbZq6BzdHQUH+L1RW5urrCaMuWj0+nEDxd6bdoTeXl5+h9nTN2+R+kzb9myZeLzkpg0aZL4Yfvjjz+yoLNB2OXKWDW333672EZHR4vt119/jbi4OCHyTMUcERQUhDlz5hiNq/qFef36dVy8eLFK+8bGxgpLFYkW+uU7Y8YMFBYWmt33559/FlZEWgf9MierIj0GU+i+yQLm6+srPthJ2G7ZssVoH7VaLSxNrVu3FvuQ4O3fvz/++OMPi+JzyFK1fft2xMTE6N3cFVmu6HoSaatWrdLvL/3CJysAjckK9vDDD8PHx0esiTh9+rTYj75EaL3BwcHiSyUtLa3MfdA5mTx5MkJDQ+Hs7CwstSTIi4qKyl1XRkYGbrrpJjRp0gSXLl0yuw+5t8aMGSP+HzRokH79prFKZBGmY9E6ab3fffddmWPRDwqyljRt2lSssVWrVvjwww9RXFyMqsbq/f777+K5pdcDvZ4tOe7HH3+MW265RTzvdHt6XZHVpTrQc0PnwfA1dvz4cTHXo0ePMj+0+vTpU8bS1LFjR7Fees6mTp1qFO4gvc46deokjksWKBJyr732Wq28r8xBz2Hv3r3FcxgRESHOr/T6NHRD0tic25PmaX8Jen9QyAa5LOl803mn15Kp614KGTlw4ABmzpyJgIAA8Rjox2dKSorRa+DcuXP4+++/9a9DKY6wqu9R8kzQsadPn272/JF4mzdvnhjT60wSc4RMJhOhJxW9p5jGC1voGKsmKipKbOmDlKAvH/pgJeFT24wfP1580JLlpCLILXLHHXcIAUjxKfRl9v3332P37t1l9qUP+scff1x8ydCHLMX+ffbZZ+KDnyyM3t7eYj/6kO/Xrx/CwsIwe/Zs8YG9bt068eW2YcMG8cVA0JcNHYeslyQ+VCqVcDX/+++/uOuuu6r8WF9//XVkZWWJL4BPP/1UzHl4eJS7Pz0+6T6feuopMUdfmIbQFx0JTXLNSueQhObVq1fFOSAxR4/zm2++EVtyBUlftPHx8eLYJAjo+CTWSeCRWCGLjqEbSiI1NVU85vT0dPG8ma5HgoQEPU+ff/65EBMUR0RIW+LKlSviNUWCcsKECVi+fLkQoiSYSLQQtI6BAweKdT399NNo1qwZDh48iFdffRUJCQlVinUi0UmWZrr9k08+KYSCJcel187IkSPxyCOPiC/lNWvWiPO+bds2DB8+HJZAQotefxSiQMck9u3bJyxFp06dEq8tLy8vISppPdLzLr0O6YfFnXfeKUQ3PS4Kj6BwCHptk2iQIPFOgvChhx4SP2boR1ZN31fmOHPmDAYPHizEFK1Po9HgrbfeKvf+qgI9HnrstHb60UBCjh4niS36AWNqaXzuuefEDxq6X9qXnjuykq1du1ZcT2Pah95r9B4kyltfee9RutDnAR2TftgaCrbVq1eL9x69PsyxaNEi8X588803q31OGCtGxzBWwIoVK0gB6P78809dSkqK7saNG7o1a9bo/Pz8dK6urrrY2Fixn4+Pj65r167Vuo+pU6eK+yiPgQMHVni9xMKFC8V+69at08/l5ubqWrVqJeb37Nkj5oqKinSBgYG6Tp066fLz8/X7btu2Tez35ptv6ufuuOMOXefOnXUFBQX6ueLiYt0tt9yia926tX6OHvvw4cMtfOQ6sSbDtRF0nPDw8Cofw93dXTdhwoQy82+99ZY49rhx48pcl5eXV2Zu9erVYv+9e/fq58aPH69zdHTUHT16tMz+dB4MXyO0T0JCgq5jx466li1b6q5du1bp2n/++ecyj1+CzoHpepKTk3XOzs66WbNm6efmzp0rzsHly5eNbj979mydk5OT7vr16xWuQbqfHTt2GM1bclzT80mvMXp93X777WXuy9xzZQq9Bm666Sb9+P777xcXut/ffvtNzP37779i3Zs3b9afG4VCoRs8eLBOq9Xqb/vFF1+I/ZYvX17mPfXVV1+VuW+6ji6Wvq/K47777tO5uLjoYmJi9HPnz58Xj8XwfR0dHS3G9Hoyhebp9VzR6/fQoUNiv++++04/J70277zzTv3rlZgxY4a4/8zMTP0cvW4NH3d13qO///672Fd6jiS6dOli9tjEDz/8oHNwcNA988wzZq9nGj/scmWsCvrFT7+wyfVEv4rp1+jGjRuF5Yogq4Gnp2ed3De5OiqzzkmJFiEhIUZWQvqlbmjBIMhylpycLFw2hrExZEkhCxS5UwiyMJEVYuzYscjOzhaWJ7qQZWPIkCGIjIzUu2jJokLWLZqzNp555pkyc4bubooXo8d18803izFZFQmyAFH27IgRI/Txk4aYZiaTxYIsWuR+JutSeHh4jdfeoUMH3HrrrfoxvQbJekbWDEPXOe1DFhjpOaILvWa1Wq1YS2WQG5meU0MsOa7h+SR3M1lw6LbSubQU6bbkTpdclpRw1K1bN2GtI2hLz4HkRv/zzz+FdZBcxGTNkyCLI1n0pNe1BLlkyUJbW+8rc9B5Ilc2WbTJwilBVljT820JhuebXm/0niR3OL0PzZ1zWqvh65XOL62N3Ka1Cb02yIJJsXASlDxBbnSygppy4cIFEepw7733YvHixbW6FsZ6YJcrY1XQhw2VK6FYD3JF0Jeq4ZcGfWGQ6GlI6MOZPtRNhQat1XQ/c/MECTr68pTcfSQk33jjDXExBwlDErWUBUwfynSOyGU2dOhQkcFGWbwNDYkVU0iskmuOXIP0GAwhMUJQjBEJdXo8VYEeL70+6EuK3Li1gaEIkCCBRaJJgkQ0fWGS2DOH6eOr6jmy5LjkWqXSPCdPnjSKLauoHA8JCsM4LoLiNMmNTYKDXJOHDh0SP6LovmiOfjQYCjoSvHSbil7XdDyKPTQVL/S6Necyr+77yhz0+MhlSy5/U+j2JBarAx2TQhxWrFghflQZ/uCTXr8VvY7oNUQYvo5qA/pMJLcquX+lJBMSd/TDUYoXNYTcsyTCKe7R8POUsS1Y0DFWBcVRmbPSGAoh+jKjD6eqfEk0BqTAd6qrV541gb7opHgwiiukrF8qzUIlWSi+hmpLUYxbQ2Iu+YSsjhSD9NJLLwmrD1lc6fGSEK1KIoE57r//fpGwQPFkUvB3TTGMQzLE8Auc1ksxey+//LLZfUlkV+ccVfW4JKwo1o1eA/TFTNYsilUjsUHlLcrjxo0bZYTknj17RBwYvddIBJAVkMQIJSLQ/ZGoo/sg0Uj3K8VwVgdry+ItT/yS8DWF4t3o/JI1sm/fvlAqleL25D0w9/qtyuuoNmN+P/roI2HdprhMeg1Q0g2t0RQpCYleM4ztwoKOaVSQW46sCZQoQB9iDQG5+Mi9QR/Shl8OplmWkiuQ5qVsXcN9pevJqkHQlzO5UiqDLCXkwqILZbzRFzwFgVsq6Cwtsmzp/mSV2LVrl7DQGQZhm7qLyTJFllc6p1WBvmRJ4NIx6cuLkkhqe+3moKQLOt9VeY7q4rj0mifxRa5FcmNKkOCoCLJimmZBd+3aVWzpRxH9iCLRRoJOcjvTlsQcWX0okYdeY+Ze19Jrl6AfWZSNXt3zU9X3lTnoNUTC0VwoguntJauZaUauObcoJeVQkswnn3xiFDpgeltLsOS1WNG+ZNHu3r27eI4oYYOSSSjpwRz0/NWFqGSsC7a9Mo0KitOiX5mzZs0yW1OOXEbkkqoOVS1bQjFGlJVpWC6C3B6UvWkIWT/I4kHWM0P32G+//SbchVJWIu1D1hIqsUBZjaYYustMy32QxYvEjSWlHSQok9ac26ii/S35IpOsFaZfJKbZoOQCotinrVu3mm3rZe6LiFzTZNGkTFByO1Vl7URNvojJ2kg/JkhQmULHJddlXR6Xzid9wRtakiiTsrLuHSQCSWQZXiRRI4m3w4cPC6udJOiovA7Fn1HpFGkfCbo9CUHKGjZ8bqjeGb2eLM22tfR9ZQ46N2TdpnNB72MJep+Znlf68UCLine truncated
"text/plain": [
"<Figure size 700x650 with 1 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"\n",
"\n",
"def validation_plot(real_quantity, pc_score, score_name, xlabel):\n",
" \"\"\"Put a PC score and the real-world quantity it should track on the same\n",
" standardized scale. If the PC *were* that quantity every dot would sit on the\n",
" dashed 45\u00b0 line; the red line is the actual fit (its slope is the correlation\n",
" r). The closer red hugs the diagonal, the stronger the validation.\"\"\"\n",
" zx = (real_quantity - real_quantity.mean()) / real_quantity.std()\n",
" zy = (pc_score - pc_score.mean()) / pc_score.std()\n",
" r = real_quantity.corr(pc_score)\n",
"\n",
" fig, ax = plt.subplots(figsize=(7, 6.5))\n",
" ax.scatter(zx, zy, s=4, alpha=0.12, color=\"steelblue\")\n",
" line = np.array([-3, 3])\n",
" ax.plot(line, line, \"k--\", lw=1.3, label=\"perfect match (y = x)\")\n",
" ax.plot(line, r * line, color=\"crimson\", lw=2.5,\n",
" label=f\"actual relationship (r = {r:.2f})\")\n",
" ax.set(xlim=(-3, 3), ylim=(-3, 3), aspect=\"equal\",\n",
" xlabel=f\"{xlabel} (standardized)\",\n",
" ylabel=f\"{score_name} score (standardized)\")\n",
" ax.set_title(f\"{score_name}: does it track the real-world quantity?\")\n",
" ax.legend(loc=\"upper left\")\n",
" plt.tight_layout()\n",
"\n",
"\n",
"park_mean = df_matrix.mean(axis=1)\n",
"validation_plot(park_mean, scores[\"PC1\"], \"PC1\", \"Park-wide mean wait\")"
]
},
{
"cell_type": "markdown",
"id": "5c88f656",
"metadata": {},
"source": [
"## PC2 \u2014 Kiddieland vs. Everything else\n",
"\n",
"PC2 is a **contrast**: some rides load positive, others negative. Reading the\n",
"extremes, the negative end is dominated by the tiny-tots rides (`kiddie_boats`,\n",
"`pony_carts`, `red_baron`...) and the positive end by the big rides in the rest of the park. So\n",
"my first guess was \"kiddie rides vs. everything else.\"\n",
"\n",
"### The category tag is a trap\n",
"\n",
"The park's own data tags each ride `Kiddie` / `Family` / `Thrill`, so I reached\n",
"for that to define the two groups \u2014 and got a muddy result. Two problems:\n",
"\n",
"1. The tag values secretly carry a **trailing space** (`\"Kiddie \"`, not\n",
" `\"Kiddie\"`), so naive equality checks silently matched *nothing*. (Lesson\n",
" learned: a `value_counts()` printout will never show you that.)\n",
"2. Even cleaned up, the tag **misfiles rides** relative to how they behave. The\n",
" big Ferris wheel (`giant_wheel`) and the large `grand_carousel` are tagged\n",
" `Family` but behave like main-midway rides; the small `s&g_carousel` is tagged\n",
" `Family` but behaves like a Kiddieland ride.\n",
"\n",
"PC2 isn't tracking the *taxonomy* \u2014 it's tracking **physical geography**: the\n",
"Kiddieland corner vs. the main midway. So let me define the two groups from\n",
"geometry the PCA never touched, and see if the contrast gets sharper.\n",
"\n",
"### Validating from the map\n",
"\n",
"Each ride has a GPS coordinate. I run k-means on **just the lat/lon** (no queue\n",
"data, no PCA) to carve the park into areas, and call the cluster with the lowest\n",
"mean PC2 loading \"Kiddieland.\" Then I build the contrast *main-midway mean \u2212\n",
"Kiddieland mean* and check it against PC2.\n",
"\n",
"(Small but important detail: **standardize the coordinates first.** On raw\n",
"degrees the clustering is unstable from one random seed to the next; standardized,\n",
"it's identical every run.)"
]
},
{
"cell_type": "code",
"execution_count": 20,
"id": "dad04b89",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"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\n"
]
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAABdEAAAKyCAYAAAA6kpdwAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjksIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvJkbTWQAAAAlwSFlzAAAPYQAAD2EBqD+naQABAABJREFUeJzsnQd4FFUXhs/uJptOEnrvvUhHpQgoAoIoCIqIgoKAgAKCAmJBBEGqFEV6F6QIIgiKIE2qVOm9JLRACgnp2Z3/+Q5O/slmN9lNNv28z7OEnZ1y587s7L3fPfc7OkVRFBIEQRAEQRAEQRAEQRAEQRAEIRn65IsEQRAEQRAEQRAEQRAEQRAEQRARXRAEQRAEQRAEQRAEQRAEQRBSQCLRBUEQBEEQBEEQBEEQBEEQBMEGIqILgiAIgiAIgiAIgiAIgiAIgg1ERBcEQRAEQRAEQRAEQRAEQRAEG4iILgiCIAiCIAiCIAiCIAiCIAg2EBFdEARBEARBEARBEARBEARBEGwgIrogCIIgCIIgCIIgCIIgCIIg2EBEdEEQBEEQBEEQBEEQBEEQBEGwgYjoQrahbNmy9OKLLzp9v9evXyedTkdLlixJdd23336by5HdySnldPb5tGjRgl/ZlV27dvG9hr+5DUe+R5mJvfcEyv7+++9TXiQrvzeo9y+//DLxPe4fLMP9JAiCIAiCIAiCIAg5BRHRhWSoIoet18GDB6XWnMz+/ftZaAoLC0vT9lFRUbx9bhRvsyOzZ8/OdmKykH0GUdSXq6srlS9fnnr06EFXr15Ntn54eDiNGTOGateuTd7e3uTh4UE1a9akESNG0O3btxPXW79+PXXt2pX35enpSVWqVKFhw4al+XkhCIIgCELOA219tC8ePHiQ1UURciiWwQ2ZGcSRVQE52TFYC+VBubIruTkwTBDSi0u69yDkWr766isqV65csuUVK1aknESZMmUoOjqaBa3sLKJDTMOPqZ+fX6rrz58/n8xmcxIRHduD7BypnZtE9IIFCyZr/DzzzDN8rxmNxiwrm5D1DBo0iBo2bEjx8fF07NgxmjdvHv3222906tQpKl68OK8DUb1Vq1Z08+ZNevXVV6lv37583/z777+0cOFC2rBhA128eJHXxWfY7s0336TSpUvzfr777jvasmUL7x/ie07hrbfeotdff53c3NyyuiiCIAiCIAiZAvpqkyZNyvazaoWMubahoaFUqFAhWrlyJb322mtSzYKQgxERXbDJCy+8QA0aNMixNZSQkMBCM4Qpd3d3yk1k1YBAZGQkeXl5ZcmxcwJ6vT7X3WsZRW6+l5o1a0ZdunTh/7/zzjtUuXJlFtaXLl1Kn3zyCT+bXnnlFbp37x5HeDRt2jTJ9l9//TVNnDgx8f26deuSNcrr169PPXv2pB9//JHeffddyikYDAZ+CYIgCIIg5BUk4ClnBLZl1LX9448/OLK7devWmVw6QRCcjdi5COliypQp1LhxYypQoABHQ0LYgeBjjRUrVlCjRo3YjsDf35+jdrdt25Zsvb///pvXgxgJ+4Jly5alWg51ehjKM336dKpQoQJHOp49e9bm1LFffvmFrRNwHPxF5Kc1IMRjnzVq1OB1ixQpQv369eMR5dRAVCmilXEe2LZo0aLUq1cvCg4OTlwHU+o+/vhj/j8i/1UriJQ8g7XT0rAeRrYBfsDV7dWpenfv3mUhr2TJklwnxYoVo5dffjlVT2IcAxYTV65coXbt2pGPjw91797doTrZuHEjtW/fnqNocWxcl7Fjx5LJZCJnERQURL179+YyoCywxoBYaQnKPGPGDKpVqxavhzpr27YtHTlyJHGdxYsX07PPPkuFCxfm8lavXp1++OGHJPtBvZ85c4Z2796dWNdqg8nW1Le1a9fydwPfEUSwI6L41q1bVusbyzt27Mj/Rxk/+uijdNUX9of6Ua8B7rH+/ftTXFxc4jqIikY0dP78+fn7+dRTT3HktD389ddfLBpDEMcsCtxb586dszr9GN/HN954g7//WuEYzwa1flAGRCoHBAQkOxYiunEPYT08I/bu3etwfUB0hh0K7gEcc8+ePYmf7dy5k8tp7VmAyBF8duDAAYePiXsKXLt2jf/+/PPPdPLkSfr000+TCeggX758LKSrWGuQd+rUif9a1rWzvzf2PuNjY2Ppww8/5HsWz4qXXnqJAgMDk61nzRNdzYdhz7Mfz9TmzZtzWfBMGzduHH9vxWddEARBEAQhcwNiHAFtNbTx8mIwBWaPNmnSxK4Z54IgZG9ERBds8vDhQ/bc07604i+AKFm3bl22fhk/fjy5uLiwGGcpwEHcxTR+jDxjXbwvVaoUC3BaLl++zBGczz//PE2dOpXFNoiLEC3tAWLKrFmz2P4A20OQswbE+86dO/OP+YQJE1i0hNCsFVRVIA5D5MYPH84X60GIa9OmDds1pMSff/7JAiW2QbkgDv70008sSiuKwusgIrVbt278/2+//ZaWL1/OL1UYTw2spwq9ENbU7bFfgPOEKIgywIYEEbERERFsI5EaiJjFeUJUhpiGfTlSJxDMIAYPHTqU14MA98UXX9DIkSPJGSCaAQIjzhcC/+TJk8nX15fvGRxPCwTDIUOG8H2HKF+UAQ05rcc/6hFREqNGjeL7B+sOGDCAvv/++8R1MHgA8a5q1aqJdQ0x1BaoA0zbQ4MR91qfPn3Y4xriqaWnNcRy1CEES9Q3xEKUA+JxWoCvNkRJ3HPw1J45cyZ/DzEAgKgJgGhoiKSIkMC5QryNiYlhEdTWwJLK9u3bubwQZCGU4zrDmgj3hbVBGjwbcFw8K1APAMeDZ3ilSpVo2rRpfI127NjBg2za+oHFCe47DERhyiSOgTJaE9ttgfPG/jGIgecQnmcYSDl9+jR/jnsJ1xz3siVYBgH/6aefJkfBQBTAdQW//vor/8W1SCsYHAMYlMnI7429z3hEw+O7gQibb775hp/1GECzF3ue/RgQatmyJS9DRD9Ee1wXyzILgiAIQl7hxo0bbLWJgCC06VJCHbRGsAdmG2NAGsElavAH2qdqsAna7MePH0+2j/Pnz/PvNfpYWA/7Uds1KiEhIRwEgn2hH4DgAMxwRgCBFjX4ZM2aNdweRPsa+3zuuee4XaDl0qVL3A9BOxDrYF30q9BfTY3MCGZBHxJtYuwfx0HQCgKn7Al4sifoShuUgrpRLUDRfkM/TG3XOxrcgPsH7X8EmKDcaKuinWfZjleDINCWxvroG+IaOBLoYhnYZplLSPuy9DDfunVrYtAOzgdtTGv6gL1BcvaAgCP0W3HvoJ5xbJQBQTfac0rp2qqBXL///nuq7WJ7gpoc+c5Yw97rnVnfpfQEKwpClqEIggWLFy+Gumv15ebmlmTdqKioJO/j4uKUmjVrKs8++2ziskuXLil6vV7p1KmTYjKZkqxvNpsT/1+mTBk+xp49exKXBQUF8TGHDRuW4nW6du0ab5svXz7extpnOC+VOnXqKMWKFVPCwsISl23bto3XQzlU9u7dy8t+/PHHJPv8/fffrS63xLJ+wKpVq5Kd5+TJk3kZymoPPXv2TFLO+/fv8/ajR49Osl5oaCgvx/4dBcfAtiNHjkyy3JE6sXb+/fr1Uzw9PZWYmBib52OL5s2b80tl+vTpfMwVK1YkuQeffvppxdvbWwkPD+dlf/31F683aNCgZPvU3oPWytumTRulfPnySZbVqFEjSTlUdu7cycfBX7UshQsX5u9EdHR04nqbN2/m9b744oskdYBlX331VZJ91q1bV6lfv76SFnr06MHfvX/++cfmeQ8ZMoSPi+uqEhERoZQrV04pW7Zs4nfW1vcI5xccHJy47OTJk3xMHFsF9yW27datW5IyXL9+XTEYDMrXX3+dZPmpU6cUFxeXxOVqPeJ4sbGxievNmzeP92vtWliiPsOOHDmSuOzGjRuKu7s7P5tUPvnkE37maJ8NeKagPJbfL1vXf9GiRfydvH37tvLbb79xPep0usTrgGvq6+urpIfevXtz3V28eDHDvjf2PuNPnDjB+xswYECSdd94441kzyX190X7rLP32f/BBx9wPR4/fjxxGe69/PnzO/T8FARBEISciNqeQhsDXL58WSldujS3j9RlKYHf2ypVqnAf6Msvv1S+/fZbpUSJEvzbjzYB9vXNN9/wC+2UihUrJum7nT59mpdXr15dmThxovLdd98pzzzzDP82r1+/PnE9tHcqVKjAfYi5c+dy2xbHwba3bt1K1m5S27ooD8qFfkKjRo0S10PbD+3S4sWLK+PGjVMWLFigjBkzRmnYsCG3JVNCbXdgXewfZfLw8OC2GfpJ2nY42oRo4/fq1Uv54YcflM6dO/O2s2fPTvEY9+7dU/z9/ZXKlStzn2v+/PnKp59+qlSrVo0/f/ToEe8P+0Kbc/ny5fxCmxlMmTJFadasGdcT2raDBw/mMqIOtP0U9fqjvl555RUu17vvvsvLhg8fnqRMb775Ji9HWwzXCes/8cQTydpla9euVWrXrs19Ehx71KhRfC64VyIjI5PVI6492pSzZs3i+wTgeuCzxo0bKzNnzuS+hZ+fH/eftO1Py77E3bt3E+tCfWG/rq6ufL1Uli1bxvdY27Zt+XPce7h+OIa27ffHH39wHwTt1GnTpvE1wD2Ha5qWfia+U/iuDB06lK/fpEmT+PuD8qlt0dSuLTh48CB/ju+PCsqDe04FdVGkSBHFx8eHy43y47rgfLTfLXu/M7aw93pb9mkz6ruE+xd9rD59+ihz5sxRRowYoXh5efEx0OcQhOyIiOhCMtQH5Pfff6/8+eefSV4QI20REhLCPzb9+/fnHzVLgVgrfFgDD2/8MFuCH3ytyGUN9Uf5nXfesfmZ+oMNYcuaOAxwfO2PLERX/PhC0MG5aV9ocOLBby8QUbGdWh4IWRktokOoNhqNSvv27fn6OIIq6kJo1JLWOoEwh3XQSMd+IbzZOh97GzetW7dWihYtmmxwRh2Line truncated
"text/plain": [
"<Figure size 1500x700 with 3 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"from sklearn.cluster import KMeans\n",
"\n",
"# category tag, with the sneaky trailing space trimmed off\n",
"category = pd.read_sql(\n",
" text(\n",
" \"\"\"\n",
" SELECT r.ride AS ride_name, trim(c.name) AS category\n",
" FROM knoebels.ride_master r\n",
" JOIN knoebels.\"LZ_attractions_io_poi\" p ON p._id = r.attractions_dot_io_id\n",
" LEFT JOIN knoebels.\"LZ_attractions_io_categories\" c ON c._id = p.category\n",
" \"\"\"\n",
" ),\n",
" engine,\n",
").set_index(\"ride_name\")[\"category\"]\n",
"\n",
"coords = pd.read_sql(\n",
" text(\n",
" \"\"\"\n",
" SELECT r.ride AS ride_name,\n",
" (p.location)[0]::float AS lat, (p.location)[1]::float AS lon\n",
" FROM knoebels.ride_master r\n",
" JOIN knoebels.\"LZ_attractions_io_poi\" p ON p._id = r.attractions_dot_io_id\n",
" WHERE p.location IS NOT NULL\n",
" \"\"\"\n",
" ),\n",
" engine,\n",
").set_index(\"ride_name\")\n",
"\n",
"annotated = loadings.join(category.rename(\"category\"))\n",
"located = annotated.join(coords, how=\"inner\")\n",
"\n",
"# k-means on STANDARDIZED coordinates only -> reproducible spatial areas\n",
"located[\"area\"] = KMeans(n_clusters=5, n_init=10, random_state=0).fit_predict(\n",
" StandardScaler().fit_transform(located[[\"lat\", \"lon\"]])\n",
")\n",
"kiddie_area = located.groupby(\"area\")[\"PC2\"].mean().idxmin()\n",
"kiddie_rides = located.index[located[\"area\"] == kiddie_area]\n",
"main_rides = located.index[located[\"area\"] != kiddie_area]\n",
"\n",
"# Does the map agree? Left: each ride at its real GPS point, colored by PC2\n",
"# loading. Right: areas from clustering lat/lon alone \u2014 no queue data, no PCA.\n",
"fig, axes = plt.subplots(1, 2, figsize=(15, 7))\n",
"sc = axes[0].scatter(located[\"lon\"], located[\"lat\"], c=located[\"PC2\"],\n",
" cmap=\"RdBu_r\", s=140, edgecolor=\"black\", lw=0.4)\n",
"axes[0].set(title=\"Each ride at its real location, colored by PC2 loading\\n\"\n",
" \"(blue = Kiddieland end, red = main-midway end)\",\n",
" xlabel=\"Longitude\", ylabel=\"Latitude\")\n",
"plt.colorbar(sc, ax=axes[0], label=\"PC2 loading\")\n",
"for ride, row in located.iterrows():\n",
" axes[0].annotate(ride, (row[\"lon\"], row[\"lat\"]), fontsize=6, alpha=0.6,\n",
" xytext=(3, 3), textcoords=\"offset points\")\n",
"for area, g in located.groupby(\"area\"):\n",
" label = f\"Area {area}\" + (\" (Kiddieland)\" if area == kiddie_area else \"\")\n",
" axes[1].scatter(g[\"lon\"], g[\"lat\"], label=label, s=140, alpha=0.85)\n",
"axes[1].set(title=\"k-means on standardized lat/lon alone\\n(no queue data, no PCA)\",\n",
" xlabel=\"Longitude\", ylabel=\"Latitude\")\n",
"axes[1].legend()\n",
"plt.tight_layout()\n",
"print(f\"Kiddieland ({len(kiddie_rides)} rides): {', '.join(sorted(kiddie_rides))}\")"
]
},
{
"cell_type": "markdown",
"id": "2bfb6369",
"metadata": {},
"source": [
"The blue (negative-PC2) rides really do huddle in one corner of the park \u2014\n",
"and clustering the coordinates on their own rediscovers that same corner without\n",
"ever touching a wait time. So PC2 is genuinely spatial. Now the quantitative\n",
"check: take *main-midway mean \u2212 Kiddieland mean* at each moment and see how\n",
"tightly it tracks the PC2 score."
]
},
{
"cell_type": "code",
"execution_count": 21,
"id": "1f8c5e7d",
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAnQAAAKACAYAAAAVauXAAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjksIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvJkbTWQAAAAlwSFlzAAAPYQAAD2EBqD+naQAA/htJREFUeJzsnQd4U9Ubxt/uvXfL3hsZgiKKOABBGe498e/eoqLiAvfeAxlOHICKgIAiIEtki2xKKdA9aNM983/eU264SdM2adP9/Z6nPTk3NzcnN+vNN52MRqMRgiAIgiAIQrPFubEXIAiCIAiCINQNEXSCIAiCIAjNHBF0giAIgiAIzRwRdIIgCIIgCM0cEXSCIAiCIAjNHBF0giAIgiAIzRwRdIIgCIIgCM0cEXSCIAiCIAjNHBF0giAIgiAIzRwRdILgIFavXg0nJyc1NjWa8tpsZe7cueoxbNmyxaHH7dChAy6++GK0JPiYbr75ZjRlzj33XPXXGl67ttBaHqdQf4igE5rUl7X25+npiW7duuHee+9FSkpKpf257dFHH0WPHj3g7e0NHx8fDBo0CDNmzEBWVpbap7y8XB13/PjxaNu2rdqnT58+ap/CwkK0dr799lu88847Nu//0ksv4eeff0ZzZM+ePXjuuedw5MiRxl6K0Mr56KOP1OdSfbxHLeHr/f7770f37t3V52TPnj3xySef1Pp4QtPGtbEXIAh6XnjhBXTs2FEJrnXr1uHjjz/G0qVL8d9//6kPJLJ582aMHTsWubm5uP7665WQI7TcvPLKK/jrr7+wYsUK5Ofn45ZbbsEZZ5yBO++8E+Hh4di4cSOeffZZrFy5En/++acSj62Bc845BwUFBXB3dzf7suB5ffDBB20WdJdffjkmTpyI5ijonn/+eWURovVKEBpT0IWGhlayoDriPWoJf8T8/fffuOGGG9Tn3/z583HXXXepH8xN3YIr2I8IOqFJcdFFF2Hw4MHq8uTJkxESEoK33noLv/zyC6655hplfZs0aRJcXFywfft2ZaHT8+KLL2LmzJnqMj8Y169fj2HDhpmuv/3229UXuibqLrjgArQGnJ2d1Yd4Q5GXl6csokLVGI1G9cPFy8urVZ0m/tDSfpwJ9fsepYdj1qxZ6vOS3HrrrejUqRO++eYbEXQtEHG5Ck2a8847T41xcXFq/PTTT5GQkKBEnqWYIxEREXj66adNgk4v5jQoCMnevXvNth89ehT79u2zaV3Hjx9XliqKFv7yfeihh1BUVGR13x9//FFZEfnFzV/mtCryMVjC+6YFLDg4WH2wU9guWrTIbJ+SkhJlaeratavah4J3+PDh+P333+2Kz6GlasmSJYiPjze5uauzXPF6irQvvvjCtL/2C59WAM5pBbv22msRFBSk1kT+/fdftR+/RLjeyMhI9aWSkZFR6T54Tm677TZER0fDw8NDWWppTSguLq5yXSdOnMCQIUPQpk0b7N+/3+o+dG9dccUV6vLIkSNN67eMVaJFmMfiOrneL7/8stKx+IOC1hK68LnGLl264NVXX1XufVtj9ZYvX66eW74e+Hq257hvvPGGek3zeeft+bqi1aU28LnhedC/xrZu3aq2DRw4sNIPraFDh1ayNPXu3Vutl8/ZPffcYwp30ODrjGEOPC4tUBRyTz75pEPeV9bgc3j66aer57Bz587q/GqvT70bknNrbk9u5/4afH/cfffdymXJ883zzteSpeteCxnhD8iHH34YYWFh6jHwsyYtLc3sNbB7926sWbPG9DrU4ghtfY/SM8FjP/DAA1bPH8Xbyy+/rOZ8nWlijri6usLNza3a95TQfBELndCkiY2NVSM/SAm/fPjBSuFTW5KTk9VIcaXnxhtvVB+0tJxUB90i559/vhKAjE/hl9lXX32lXLiW8IOebl9+yfBDlrF/7777rvrgp4UxMDBQ7ccP+bPOOgsxMTF44okn1Af2Dz/8oL7cFixYYBKh/LLhcWi9pPgwGAzK1bxt2zZceOGFNp+Dp556CtnZ2eoL4O2331bbfH19q9yfj0+7z//9739qG78w9fCLjkKTrlntHFJoHj58WJ0Dijk+zs8++0yNdAVpX7SJiYnq2BQEPD7FOgUexQotOno3lEZ6erp6zJmZmep5s1yPBoUEn6f33ntPiQnGERFtJIcOHVKvKQrKm266CbNnz1ZClIKJooVwHSNGjFDruuOOO9CuXTts2LABU6dORVJSkk2xThSdtDTz9rQWUyjYc1y+dhgTet1116kv5e+++06d98WLF2PcuHGwBwotvv4YosBjkrVr1ypL0c6dO9Vry9/fX4lKrkd73rXXIX9Y0MJN0c3HxfAIhkPwtU3RoEHxTkF49dVXqx8z/NFV1/eVNXbt2oVRo0YpMcX1lZaWKkt8VfdnC3w8fOxcO380UMjxcVJs8QeMpaXxvvvuUz9oeL/cl88drWTff/+9up5z7sP3Gt+DpKr1VfUe5R8/D3hM/rDVC7Z58+ap9x5fH9Z4//331fvxmWeeqfU5EZowRkFoAsyZM4cKwPjHH38Y09LSjMeOHTN+9913xpCQEKOXl5fx+PHjar+goCBj//7963RfF1xwgdHf39944sQJs+0jRoxQa6iJd955R+33ww8/mLbl5eUZu3TporavWrVKbSsuLjaGh4cb+/TpYywoKDDtu3jxYrXfM888Y9p2/vnnG/v27WssLCw0bSsvLzcOGzbM2LVrV9M2PvZx48bZ/Zi5Jv3aCI/Tvn17m4/h4+NjvOmmmyptf/bZZ9Wxr7nmmkrX5efnV9o2b948tf9ff/1l2nbjjTcanZ2djZs3b660P8+D/jXCfZKSkoy9e/c2durUyXjkyJEa1/7jjz9WevwaPAeW60lNTTV6eHgYH3nkEdO26dOnq3Nw4MABs9s/8cQTRhcXF+PRo0erXYN2P8uWLTPbbs9xLc8nX2N8fZ133nmV7svac2UJXwNDhgwxzS+99FL1x/v97bff1LZt27apdf/yyy+mc+Pu7m4cNWqUsayszHTbDz74QO03e/bsSu+pTz75pNJ98zr+2fu+qoqJEycaPT09jfHx8aZte/bsUY9F/76Oi4tTc76eLOF2vp6re/1u3LhR7ffll1+atmmvTX62aK9X8tBDD6n7z8rKMm3j61b/uGvzHl2+fLnaV3uONPr162f12OTrr782Ojk5Ge+8806r1wvNH3G5Ck0K/uLnL2y6nvirmL9Gf/rpJ2W5IrQa+Pn51fr4tB798ccfKnlCs45p0NVRk3WOMEkjKirKzErIX+p6Cwah5Sw1NVW5bPSxMbSk0AJFdwqhhYlWiCuvvBI5OTnK8sQ/WjZGjx6NgwcPmly0XDOtW9zW1GDiiSX6+DDGi/FxMUmF0KpIaAFi9uwll1xiip/UY5m4QosFLVp0P9O61L59+zqvvVevXjj77LNNc74GaT2jNUPvOuc+tMBozxH/+JotKytTa6kJupH5nOqx57j680l3My04vK12Lu1Fuy3d6ZrLkglHp512mrLWEY58DjQ3Ot8/tA7SRUxrngYtjrToaa9rDbpkaaF11PvKGjxPdGXTok0LpwatsJbn2x7055uvN74n6Q7n+9DaOeda9a9Xnl+ujW5TR8LXBi2YjIXTYPIE3ei0glrC8BKGOkyYMAEffvihQ9ciNB3E5So0Kfhhw3IljPWgK4JfqvovDX5hUPTUBrooGF9HtxrdRLWFH878ULcUGlyr5X7WthMKOn55au4+Cslp06apP2tQGFLUMguYH8o8R3SZjRkzRmWw9evXD40NxYolFKt0zdE1yMegh2KEMMaIQp2Pxxb4ePn64JcU3biOQC8CNCiwKJo0KKL5hUmxZw3Lx2frObLnuHStsuzOjh07zGLLqsvWpqDQx3ERxmnSjU3BQdcks7/5I4r3xW380aAXdBS8vE11r2sej7GHluKFr1trLvPavq+swcdHly1d/pbw9hSLtYHHZIjDnDlz1I8q/Q8+7fVb3euIryGifx05An4m0q1K96+WZEJxxx+OWryo5WcfRTjjHvWfp0LLQgSd0KRgHJU1K41eCPHLjB9OtnxJaDCWizFytI41tTpMWuA76+pVZU3gF50WD8a4Qmb9sjTL559/ruJr+JgY49aYWMvWpNWRMUhTpkxRVh9aXPl4KURtSSSwxqWXXqoSFhhPpgV/1xV9HJIe/Rc418uYvccee8zqvhTZtTlHth6XwoqxbnwN8IuZ1izGqlFssLxFVRw7dqySkFy1apWKA+N7jSKAVkCKESYi8P4o6ngfFI28Xy2GszY0tSzeqsQvha8ljHfj+aU18swzz0RAQIC6Pb0H1l6/tryOHAU/z15//XVl3WZcJl8DTLrhGi3RkpD4mhFaLiLohGYF3XK0JjBRgB9itrBp0yb1hcQvLyYa0LpTF+jio3uDH9L6LwfLLEvNFcjtWraufl/telo1CL+cbSmjQksJXVj8Y8Ybv+AZBG6voLO3Bp+9+9MqwdIwtNDpg7At3cW0TNHyynNqC/ySpcDlMfnlxSQSR6/dGky64Pl2dKkbW4/L1zzFF12LdGNqUHBUB62YllnQ/fv3VyN/FPFHFEUbBZ3mduZIMUerDxN5+Bqz9rrWXruEP7KYjV7b82Pr+8oafA1ROFoLRbC8vWY1s8zIteYWZVIOk2TefPNNs9ABy9vagz2vxer2pUV7wIAB6jliwgaTSZj0YA0+f/UhKoWmhdhehWYF47T4K/ORRx7BgQMHKl1PlxFdUhp0y9Eqx3R/uquqsxbYWraEMUbMytSXi6Dbg9mbeiggafGg9UzvHvvtt99M6yLch9YSllhgVqMleneZZbkPWrwobuwp7aDBTFprbqPq9rfni0yzVlh+kVhmg9IFxNinX3/91WpbL2tfRHRN06LJTFC6nWxZO6nLFzGtjfwxQUFlCY9L12V9Hpfnk1/weksSMylr6t5BEUiRpf/TRI0m3vijh1Y7TdAxA5zxZyydou2Line truncated
"text/plain": [
"<Figure size 700x650 with 1 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"contrast = df_matrix[main_rides].mean(axis=1) - df_matrix[kiddie_rides].mean(axis=1)\n",
"validation_plot(contrast, scores[\"PC2\"], \"PC2\", \"Main-midway \u2212 Kiddieland wait\")"
]
},
{
"cell_type": "markdown",
"id": "efbf0dda",
"metadata": {},
"source": [
"r \u2248 0.78, built from coordinates the model never saw, and `giant_wheel` /\n",
"`grand_carousel` correctly land on the main-midway side while `s&g_carousel`\n",
"lands in Kiddieland \u2014 exactly the rides the category tag got wrong. PC2 is\n",
"**where in the park you are**."
]
},
{
"cell_type": "markdown",
"id": "d896a1de",
"metadata": {},
"source": [
"## PC3 \u2014 morning school groups vs. the evening crowd\n",
"\n",
"PC3 was the one I got wrong at first, so this is the honest version.\n",
"\n",
"My initial hypothesis was *ride throughput*: \"people-eater\" walk-on rides vs.\n",
"low-capacity rides that build a line. It's a tidy story, and it's **wrong** here:\n",
"\n",
"- The `capacity` and `duration` columns I'd have needed are **completely empty**\n",
" in this dataset, so I couldn't even compute throughput.\n",
"- When I split rides by their PC3 sign and compared actual wait behavior, the two\n",
" groups were a near-tie (mean wait 9.8 vs. 8.8 min). PC3 simply doesn't separate\n",
" \"builds a line\" from \"walk-on.\"\n",
"\n",
"What it *does* separate is **time of day**. The PC3 *score* runs strongly\n",
"positive at midday and strongly negative after dark. So the rides loading\n",
"positive are busy in the daytime, and the ones loading negative are busy in the\n",
"evening \u2014 and that maps onto a real thing about how Knoebels operates:\n",
"\n",
"> Knoebels runs discounted group rates for schools and daycares, who arrive at\n",
"> open and leave by mid-afternoon. Those big groups of little kids hammer the\n",
"> daytime family flat-rides (tea cups, tilt-a-whirl, paratrooper). After they\n",
"> clear out, the evening belongs to teens and date-night couples on the coasters,\n",
"> bumper cars, and the Ferris wheel at dusk.\n",
"\n",
"**The test:** de-trend each group's wait by the park-wide average for that hour\n",
"(so we're looking at *shape*, not PC1's overall nightly ramp), and watch the two\n",
"groups pull apart across the day."
]
},
{
"cell_type": "code",
"execution_count": 22,
"id": "a1723720",
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAABW0AAAHqCAYAAAB/bWzAAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjksIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvJkbTWQAAAAlwSFlzAAAPYQAAD2EBqD+naQAA8RVJREFUeJzs3Qd4k1UXB/B/d+mmQCelUHYZZe89BEGGAiKyUQQUB6gI+ikgKoqiKChDBQRElojgQGXvvfcqlNEBFNrS0p3vORdT05K2aUmbpP3/nieQvHmb3NykzX3Pe+65VhqNRgMiIiIiIiIiIiIiMgvWpm4AEREREREREREREf2HQVsiIiIiIiIiIiIiM8KgLREREREREREREZEZYdCWiIiIiIiIiIiIyIwwaEtERERERERERERkRhi0JSIiIiIiIiIiIjIjDNoSERERERERERERmREGbYmIiIiIiIiIiIjMCIO2RERERERERERERGaEQVsiMjkrKyuMHj26QJ9j4cKF6nkOHDgAS5aeno6aNWviww8/LNTn1fbf5cuXUdi2bNminlv+J/Mn79WkSZMK9DnmzJmDcuXKISkpqUCfh4iIqCCVL18eQ4YMKZadXBjj/6ykr6XPcyPjXWmfjH8tQUGNlaW/XFxcUBQ/C0SWgkFbojwErLQXR0dHVKlSRX25REZGPrS/bHvjjTdQrVo1ODk5wdnZGfXr18cHH3yAu3fvZuz37bffonXr1vD29oaDgwMqVKiAoUOHmiQwRpbhp59+wtWrV4vkwOabb76xmMGx1syZM+Hu7o6UlBRTN6VYkYOI5ORkzJ0719RNISIiAKGhoWpsIuNjGfvKJTg4GC+99BKOHTvGPjLQqVOn1IlPHgvQo1q6dClmzJjBjiSycLambgCRJXn//fdVYDUxMRE7duzA7Nmz8ccff+DEiRNqcCr279+PLl264N69exgwYIAK1grJ8Pz444+xbds2/P3332rb4cOH1eN1794dJUuWVANeCeT+9ttvOHr0KPz8/Ez6esn8fPrpp3jmmWdUoLAoBm1Lly79UMZJq1atcP/+fdjb28Pc/P7773jsscdgZ2dn6qaYDXmvbG0LdnghJ84GDx6Mzz//HC+//LI6mUZERKYh49a+ffuqv/39+/dHSEgIrK2tcebMGaxevVqNl2WMGxgYyLcoi7Nnz6q+0g3aTp48GW3atDEoI5TyRo6zZNZaUaNvrCxBWzlGfe2110zaNiJ6NAzaEuXB448/jgYNGqjrzz//PEqVKqWCBr/++iv69eunsmiffPJJ2NjYqICsZNrqkintMljQDVJl1bNnT/UcixYtwvjx4836/ZFBj2S7SQCF8kYC/zKw0h2o50Y+UxLMnz59eq77agf7lpa5qo/0kTl+xhISErB161Z1MEr/Kaz36umnn8a0adOwefNmtGvXjm8BEZEJXLx4UZ1MloDsxo0b4evrm+n+Tz75RI13cxvvxMfHq5lpxY3MtCtONBqNGgOXKFHCJM9fVE+ym+tY2RKZ+jNKlBXLIxA9Am2gQLIHhEzVvX79ugrkZg3YCimD8L///S/Hx9SeVdcto5Cdf/75By1atICHh4eqN1S1alW8/fbbmfaRLx2ZZiXT1eTLXAbTTz31lBpk6w6UX3/9dQQEBKjBozzOZ599pr609NUb+vHHH1GjRg217/r169V98rqHDRuWUepB7p8/fz7yQh5XnlvaKRnKkpWsJYEZef5ffvnloZ+TM8ly3+7du3N9DqmBOXbsWJQpU0YdHEiQ/ebNmw/tJwcY2tcoGc8yvS/re5JdHTIJmMola52pZcuWqfff399fZWbHxsYiL9asWaMCvXI2vSCdPHlSfbZlsFK2bFlV1kNfVoKcrOjatavqH+mnihUrYsqUKUhLS8vYZ+LEiWqArK+PX3jhBfXZlc+o9KU8rwRBtWVItH2or06X3Ce1fWXKpZQYkf6sVKkSVq1ape6Xx2ncuLF6DfKZ2rBhw0PP/6ifWTk4lc+TnMwpTEuWLFG/H/LaPD091cGylMzQkt9R+XsgQeWs5OSSj49Ppvfozz//RMuWLdXvg6urq3pP5b3QV9NM+kxOLMl1+R2SMjC6j6Wvpq1cl20XLlxQjyPvuWSKSymYrG2ULJFXXnlFZVxLW2QWgjynvjq50gfy+uVzSEREpiEnz2QcuWDBgocCtkKyb+Xvuowxs36nyFhUZqfJ33vJ0DV0TJpTvdHsvoMk61dO9rm5uamki1dffVWNP3Ly1VdfqUQM3fGfnDiXx5OxpJZ8D8preOuttzK2SZubNWumnku+r+U7SztGyW4sKa+nT58+6nrbtm0zxkN5qVOq7Zv8ZOnmpa/k/ZaxopeXl3qfpBSGvpPY0o4nnngCf/31l0pKkb7IqbSRjDklACnlp7Ij74e8L/L+aN26dUv9nLRX97MyatQoNe7JqaatPJ5sl7GJjFFkJk92x2HSN71791bjDzlekde0du1a5KZevXrq+EtXrVq1VH/rlg9Zvny52nb69Gl1+8qVK3jxxRfV74H0nbw++YxkLZ+Rdaws42SZDSY/r/0cPUrm9qVLl9CpUyc1VpRxv8w+1faz/C+P3aNHj4d+Tj430q8jRoww+FhHxvfacbn2ODNrEouMveXzKX9H2rdvjz179uj9LBuyRkdeP6NEhY2ZtkSPQBv4lC9QIV/a8odevszz4vbt22rAFxYWpr4EhXwB5USCKvIFU7t2bfUz8uUmQZGdO3dm7COPKftIcEkCOzLoiouLU8FemS4jQTb5opXAiARFn3vuOdSpU0d9ab355psqWPLFF19ket5NmzZhxYoVKjAkgRX5opMavk2aNMkI6kowRwJB8ngSmDRkWo4E2WSgIgN7eS0SNO3cuTP27dunvrxl8CEDeAnsSqBVl2yT19K0adNcn0emUkspCgkmyhe21HqSNstz637Ry9S0Dh06qMGeTF2TgaiUvpD+ze9ZegloStBVAl0S7MvrdP9du3apvijILIGIiAh1oJCamqoyvWVwNm/ePL1nm2XgI4MlOXCR/+Wz8d5776n3XMo4iIEDB6rPp/Svbh1eydCWg5devXqpQa+8D/LeyOO88847ah8Jpubkzp076vMtn20ZwMp7JNfl8yCfuZEjR+LZZ59VbZHfSQlsykGVMMZnVkqjyEFYbu00JsnWf/fdd9XBlGT7SzBcDmwkkC+DWDnYkCmqX3/9tRqsaw/+hARI161bpw5M5GBHLF68WB2cyEBcsqFkH+lHORkkj6c7wJe/J7KfBMPlYFQC4XLwKr978nuSG2mzlIOZOnUqDh06hO+++04d7Mnzaknb5O+LfG7k/ZG/CxJEzukgSPdvHhERFX5pBDlpKt8NeSHjDPlOke8b+U6Rk695HZPmhXwHyXeafAdJgEcCfjKOkJlt2ZETmnLSWkqiyXhDbN++XQUH5X8t+b6Usmi6J9W//PJL9VokGC1jHjlxL9/J0l/Zfa/Jz8s4WNomSRjVq1dX27X/FxZD+krGChJUk9cogXkZX0hwUfpLEh10yThaThpL4G748OEqAKmPJDZ89NFHKmAm+2VHxjoyHpbkDukvIe+RjOmio6NViQlpm5D3Sd7H7MhnToKN8vMybpS+lgQRGRvpO/Zq3ry5Sr7QjpFlzCIns3/++eeHjk90SRtkXQotaac8nvazJMdz2vbKmFT7nsuxh4z/ZXwriRRy7CJ9L8dF8jq15fmykrF0TEwMrl27lvF7k98FxWT8J8dkMi6TkzQSSJXjKPkdljG+9LuUBJT75HVJQFtLPhcyrpb7cyPvgZRTkc+RjNflcyfHCXJ8rD3Wlj6TvpSA7bhx49QxkXxepD+0CRv5YehnlMgkNESUqwULFsipRM2GDRs0N2/e1Fy9elWzbNkyTalSpTQlSpTQXLt2Te1XsmRJTUhISJ571MHBQT2+XOQxv/rqq1x/5osvvlD7S3uyM3/+fLXP559//tB96enp6v81a9aofT744INM9/fu3VtjZWWluXDhQsY22c/a2lpz8uTJTPs+99xzGl9fX82tW7cybX/mmWc07u7umoSEhBxfi/a1HzhwIGPblStXNI6Ojponn3wyY9uECRNUX929ezdjW1RUlMbW1lYzceJEg97DDh06ZLx2MWbMGI2NjU3GY8rj2dvbax577DFNWlpaxn6zZs1SPy99qhUYGKgZPHjwQ8/VunVrddHavHmz+tmgoKBc+yInZcuW1fTq1cugfeX59bUtN6+99ppq6969ezO2SZ/I+yjbQ0NDM7brey0jRozQODk5aRITEzO2NW3aVNO4ceNM+61evVo9nvSNVo0aNTL1W9b+091X9pNtS5cuzdh25syZjM/onj17Mrb/9ddfart8Boz1mRXlypXL9XNnTJcvX1af1Q8//DDT9uPHj6vfAe12+Xz7+/s/9FlZsWKF6odt27ap23FxcRoPDw/N8OHDM+0XERGh+kB3u3yW5Gfff//9TPvWrVtXU79+/UzbZD/dfpHrsm3YsGGZ9pPfbfl7p3Xw4EG1n3wGdQ0ZMuShx9R64YUX1N9gIiIqfDExMervc8+ePR+6786dO2qMqr3ofq9qv1PGjx+f6WcMHZPKWCTr93pu30Hdu3fPtN+LL76oth89ejTb1yfjQDc3N824ceMyvl/le6tPnz7q+1i+R4WMs2XsIa9ZK+s4Ijk5WVOzZk1Nu3btMm3POpZcuXLlQ2OevJDnOX36dKbxu6Hy0lf6xkmdOnVSY92sr09+dv369Q/tL9tfeukldf31119Xfbhw4UKD2io/5+3tnXF77NixmlatWmm8vLw0s2fPVttu376tPjdffvllxn7S19KmrJ+5adOmZWxLTU3VtGzZ8qHPWPv27TW1atXKNMaVz0SzZs00lStXzrG92vf11KlT6vbatWvVMY30dd++fTP2q127dqZjH339vHv3bvVYixYtynGs3LVr10yvNT+0v6svv/xyptcsjy3HS9rLine truncated
"text/plain": [
"<Figure size 1400x500 with 2 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"daytime_rides = loadings.index[loadings[\"PC3\"] > 0.1]\n",
"evening_rides = loadings.index[loadings[\"PC3\"] < -0.1]\n",
"\n",
"hour = df_matrix.index.hour\n",
"park_by_hour = df_matrix.mean(axis=1).groupby(hour).mean()\n",
"\n",
"fig, axes = plt.subplots(1, 2, figsize=(14, 5))\n",
"\n",
"# left: the axis itself, by hour\n",
"scores.groupby(scores.index.hour)[\"PC3\"].mean().plot(ax=axes[0], marker=\"o\")\n",
"axes[0].axhline(0, color=\"grey\", lw=0.5)\n",
"axes[0].set(title=\"PC3 score by hour (+ daytime / \u2212 evening)\",\n",
" xlabel=\"Hour (local)\", ylabel=\"Mean PC3 score\")\n",
"\n",
"# right: each group's wait as a share of the park-wide level, de-trended\n",
"for name, grp, color in [(\"daytime rides (PC3 +)\", daytime_rides, \"tab:orange\"),\n",
" (\"evening rides (PC3 \u2212)\", evening_rides, \"tab:blue\")]:\n",
" (df_matrix[grp].mean(axis=1).groupby(hour).mean() / park_by_hour).plot(\n",
" ax=axes[1], marker=\"o\", label=name, color=color)\n",
"axes[1].axhline(1.0, color=\"grey\", lw=0.5)\n",
"axes[1].set(title=\"Group wait \u00f7 park-wide wait, by hour\",\n",
" xlabel=\"Hour (local)\", ylabel=\"Share of park-wide mean wait\")\n",
"axes[1].legend()\n",
"plt.tight_layout()"
]
},
{
"cell_type": "markdown",
"id": "429f01a2",
"metadata": {},
"source": [
"The daytime group fades through the afternoon while the evening group\n",
"climbs \u2014 and the school-group operating detail is part of the mechanism. PC3 is **what\n",
"time of day it is**.\n",
"\n",
"(There's a mild bit of circularity worth flagging: I defined the two groups using\n",
"the PC3 loadings, so this plot confirms *what the axis means* more than it's a\n",
"fully independent test. The genuinely external evidence is the operational fact\n",
"about group rates.)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": "### What PC3 adds that PC2 doesn't\n\nIt's fair to ask what's left for PC3 to do. PC2's *score* also drifts with the clock \u2014 kiddie-heavy in the morning, midway-heavy by night \u2014 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.\n\nThe mechanics: each component assigns every ride a **single number** (its loading), and the axes are built to be independent \u2014 here PC2 and PC3 correlate **0.00**. PC2 spends its one-number-per-ride on a single question \u2014 *kiddie corner or main midway?* \u2014 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."
},
{
"cell_type": "code",
"metadata": {},
"source": "# PC2 calls these two rides almost identical (both lean \"main midway\").\n# PC3 sends them to opposite ends.\nloadings.loc[[\"tilt-a-whirl\", \"giant_wheel\"], [\"PC2\", \"PC3\"]].round(3)",
"execution_count": null,
"outputs": []
},
{
"cell_type": "markdown",
"metadata": {},
"source": "PC2 puts them practically side by side (+0.16 and +0.20); PC3 sends them to opposite poles (+0.15 vs \u22120.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 \u2014 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.\n\nIt 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:\n\n- **PC1** \u2014 *how many* people are here? (overall busyness)\n- **PC2** \u2014 *which half* of the park are they in? (kiddie corner \u2194 main midway)\n- **PC3** \u2014 among the midway rides, *which kind* is peaking? (daytime family flat-rides \u2194 evening marquee rides)\n\nThe gut-check for what \"independent axes\" really buys you: two moments can share an identical PC1 **and** PC2 \u2014 same crowd size, same kiddie/midway balance \u2014 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."
},
{
"cell_type": "markdown",
"id": "ba034336",
"metadata": {},
"source": [
"## The map of the park, drawn from wait times alone\n",
"\n",
"Here's the payoff. Plot every ride by its PC2 and PC3 loadings \u2014 its *behavioral*\n",
"position \u2014 and color it by the park's category tag. Remember: the **position**\n",
"comes purely from how each ride's queue moves over two seasons. The model was\n",
"never told what any ride is or where it sits.\n",
"\n",
"And yet the picture is the park. Kiddieland collapses into a tight cluster on the\n",
"left. The evening thrill rides sink to the bottom-right. The daytime family\n",
"flat-rides float to the top. The two blue dots stranded inside Kiddieland are the\n",
"`Family`-tagged rides the queue data insists are really kiddie rides."
]
},
{
"cell_type": "code",
"execution_count": 23,
"id": "3a436334",
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAABD0AAANHCAYAAAA16iZRAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjksIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvJkbTWQAAAAlwSFlzAAAPYQAAD2EBqD+naQABAABJREFUeJzs3Qd4VGX2x/GTSQ9JCEhHipQFRSmCBRu2tVdUFP3bdXftvfeK61pXLGtdK+La3XXtiI1VQFREwIbSBaWEENLn//ze4cZJmElmwkzKzffzPAOZO29m7ty5M5l77nnPSQkGg0EDAAAAAADwmUBTrwAAAAAAAEAyEPQAAAAAAAC+RNADAAAAAAD4EkEPAAAAAADgSwQ9AAAAAACALxH0AAAAAAAAvkTQAwAAAAAA+BJBDwAAAAAA4EsEPQAAAAAAgC8R9AAAtGq77rqrbbnllgm9z969e9sJJ5xgyfTPf/7TUlJS7KeffmoW61Ob1u3aa69t1McEAACojaAHACCpB+XTpk2rsXz16tW27bbbWlZWlr3xxhtsfcREQRvtT94lPz/fhgwZYrfffruVlpZuMP6LL76w//u//7MePXpYZmamtW/f3vbcc0977LHHrLKy0o357bff7G9/+5vtsssu1rFjRysoKLDtt9/eJk6c2CSvSvjzCwQC1q1bN9trr73s/fff32CsnoOei4J2em56jgpunXjiiTXec7NmzbIjjjjC+vTpYzk5OdahQwf3fF977bVGfnYAADSNtCZ6XABAK1RYWOgO4r766it76aWXbJ999mnqVWqxjj32WDvqqKPcwW5ztG7dOktLS+zXDD3Xhx9+2P28atUqe+GFF+zCCy+0qVOn2rPPPls9TmP+8pe/WOfOnd126t+/v61Zs8beffddO/nkk23JkiV2+eWX25QpU+yKK66w/fbbz6688kq3vrpPbddvvvnGrrvuOmtsf/zjH+24446zYDBo8+bNs/vuu8923313+89//mP77rtv9bYdPXq0CxoqgKHnosCHsn6ee+45e/zxx23+/Pm26aab2s8//+ye+/HHH++CKMXFxe45HnTQQfaPf/zD/vSnPzX6cwQAoFEFAQBIgsceeyyoPzNTp0511wsLC4Pbb799MCMjI/jvf/+72WzzUaNGBQcNGpTQ++zVq1fw+OOPDyZDUVFRs1qfxqL1b9OmTY1llZWVwREjRrj9bNGiRW7ZlClTgqmpqcGddtrJ7XO1aX/Uvik//vhj8Keffqpxe1VVVXD33XcPZmZmxrStNf7RRx91/0fywQcfBGfPnh3Tc9TzOOOMM2os++qrr9zyvfbaq3qZxmjZnXfeucF9VFRUBP/2t78FFyxYEPVxNGbIkCHBAQMGxLReAAC0ZExvAQAkXVFRkcvq+Pzzz91Z5v3333+DqQu5ubm2aNEiO+SQQ9zPmm6gs/jeVATP2rVr7YILLqietjBgwAC77bbb3Jnx2p566ikbPny4ZWdnuzPhOoO/YMGCiOs4ffp022GHHdzYzTbbzB544IENxmgaxTXXXGP9+vVzj611uPjiiyNOrwhXXl7usgaUcaBpPZtssonttNNO9vbbb8c0RWjy5Ml2+umnW6dOndzZ+2g1PbQNbrzxRjdGUxl22203N70hEmVKnHvuudXbUc/pr3/9q1VVVdUYpwwKbcO8vDw3pWSrrbayu+++2+Kt6aGftez77793r7emkrRt29ZNx1D2QUNoCoimd4i3HbSd9ThPP/20W+faRowYUV3fRK9zr169Nlhv7YN6TX/88cd61+Hjjz+2k046ye2Ttc2YMcMOOOAAu/rqq62htL01JUVZH7Jw4UKXoaGMEL1+taWmprr3jbefRKIxet21D9Rn6dKl7jXS/Wk/6dq1qx188ME19jtNq9HzfOutt2zo0KFuH99iiy3sxRdfrHFfK1ascOum56T3uPYnZa98+eWXGzxuSUmJ22f+8Ic/uPvT4yq75Ycffqgeo331rrvuskGDBrkxyuz585//bCtXrqz3eQEAWg+mtwAAkkpBCh3YaArC888/7w6OIlFwY++997btttvOBTHeeecdV6+hb9++dtppp1Uf1Cstf9KkSW6agg6w3nzzTbvoootcwOTOO++svr+bbrrJrrrqKhszZoydcsoptnz5crvnnnvcdAAdjOqg26ODJE1x0NixY8e6KQJ6zIyMDHdA6x1g6bE/+ugjNyVg8803t5kzZ7rH/Pbbb+3ll1+Oug108DZu3Di3Hqpnomk+qrugIJAOXuujgIeCQDp41vaMRrcr6KHnoovuX9OJysrKaoxTkGHUqFFum+kgsWfPnvbJJ5/YZZdd5qZ+6EBSFJTR9thjjz1cQERmz57tDvTPOeccawhtYwUbtD20fpqKomCOd//x8g6CFUjS89IUFr3Gek4NpQN9UbChPgpe3XrrrS741a5dO7fPyXfffecCfVoPBSkaSvumLgpKyX//+1+rqKhw03biof1G02JUU+fVV19193PkkUfW+3uHHXaYC5ydddZZLrixbNkyt19o+oyue/R8dX+aVqSpNKo3oloimoLj7eMKIul9ouXaB3755Re3bbQvajqRpt94nwX6nNBrqUCl9jVN0dHjfv311+4zQbTvKvinoMzZZ5/tAkPjx49372/to+np6XFtIwCATzV1qgkAwN/TWzS1Ij09Pfjyyy/XOXVBY6+//voay4cNGxYcPnx49XXdh8bdeOONNcYdfvjhwZSUlOD333/vrmvKgqY43HTTTTXGzZw5M5iWllZjuaa36D5vv/326mWlpaXBoUOHBjt16hQsKytzy5588slgIBAIfvjhhzXu84EHHnC///HHH0edTqKpBPvvv3+wodtQUzU0JSHSbfPmzXPXly1b5qYO6XHCp1pcfvnlblz4+txwww1uqsi3335b4z4vvfRSt93mz5/vrp9zzjnB/Pz8DR47FnrMa665pvq6ftayk046qca4Qw89NLjJJpvEPL1l+fLl7qLX+uabb3av++DBg92YL7/80j2G1ruhfvvtN/e677zzznH93iWXXOIee/z48cGFCxe6fWCzzTYLLl68OOb70O+ffPLJ7vnp9fz000+De+yxR43987zzznPXZ8yYEdf6/fnPf3a/p4v2Y71nVqxYUefvrFy50o3XdJm66Llq3AsvvFC9bPXq1cGuXbu697CnpKTETUkKp/1XU4nC3/uaLqT7u+OOOzZ4LG/f1vtQY55++ukat7/xxhsRlwMAWi+mtwAAkkpnc5V6rnT6+ugscbidd965xhSD119/3aXm66xuOE0t0DGjzl6L0uqVmaGsgl9//bX60qVLFzfFRJki4VTAUmeNPcrw0HWd1da0F/nXv/7lsjsGDhxY4z5VZFJq32c4ZZXobLnOhjfEqaee6p53XZQZo4wOnZHXFA1PpCkQei7atspMCH8u6m6is+wffPBB9XorQ6C+aTjxiPQaq4uKsl/qo3VRxosuynxQAc+RI0e6orji3UekaS2x0D5zzDHHuGkfygqKxy233OIyebT91QFG02O03TQtIx6PPPKIe37KflHWkzIWzj///OrXsaHPUb+v9VGRU2Ve6XWunQFUm6Z66b2g7jH1TRlRlsahhx5afV1TV1SQVVkXXuaMpsdoSpLo8fW6a5qLpqgp68ejKXDKstG2rM3bt7UPa3qUskjC92FNxdJ91vV+BAC0LkxvAQAkldLXddCmVP8PP/zQHeBEosCIDvbC6aA8/GBLnSh0cFX7gE/BCO92UXBBQRAFOCKpnfau+2zTpk2NZaolIKpdoINY3aemdtReR48CJNFcf/31rg6C7nPLLbd020LTEwYPHmyx0FSA+njPvfZz1vpqO4bTc1EHnfqei6bVaKqPDpK7d+/upsookLQxXXdqTzvx1k2vsw6U66J9xGu1qgNobZfw2hXe72sqREPoIFvTMZ544gnXDjdean+rg3HV3bjjjjuqp2HEQ/vJmWee6Q7utZ+rXkX4vtnQ56hgnS6iYIReywMPPNA+/fTTGkGycNrGmnakoKLqZeh9oGkn+n0FEMMpCFX7fsLfQxqvoJLqwagjjaaihNfr0fSk8ClL+pyoq/uP9mFN1VFwKN73IwCgdSHoAQBIKhU0VIaG6kLorKzOXEfK+qgvkyEeOrjSAZgyPyLdr84EN+Q+VYBRB7OR1JXJohoTOpB75ZVXXLFH1bFQLRAVS1V2QH10xj2R9Fz0WqgORSTewaoOKL/44gtXN0XbUhfVatBBrzIGGiLa6xypEG2k31U2SjQ68NaBsmqtxEsFUHUwroyNeOtleIU3VQBV9TZ23HFHl4UybNiw6kKrsVIQp67n6AUu9BxV06ahDj/8cJfNpHo00QKRXoaIgiOqxaH9QDVLVI/lvffec88vHjfffLP7fdXJueGGG1xxYWV+6DFqF9Ctj8Zr/1TB2kiiBfQAAK0PQQ8AQNKpeKcOmtS1RQfbyvhoyEGJOm1oGofOcodne8yZM6f6dtEZdh1EKxPAO4Cvy+LFi93UifAz6joYFK9Yo+5TXSYUvIl2ZrwuOsBTwUVd1M1GgRAVOI0l6BEL77nrDHifPn2ql6uAa+2pCXouWoe6Dq49mt6gg15ddKCp7A9l7+jg1Suu2VyoY42mG+mAXF16YplSJffee697LXTwfckll8T9uAp0KAPmf//7n/373/92XYCUSeEV3dWUi0RR1o2CP+pM1JDgjEdFTUXZEvXR/qJsD120fynYoiLDWgePuvLoPRf+3qj9HlIhY3UU0hSecJpOFF40Vo+nDBR1PYpWjFRj9FmgAFOig4IAAH+hpgcAoFEoWDBhwgR3cKTpEbHUcKhNHUmUEq8ODeGUNaGDLR0Qilpb6sBQZ+9rZxDoumoJ1D5oDe+woVoHuq7AjHfAqoNadTt56KGHIh5A1tVVpfbjKdNEAYP6Wt3GQwEMHSCqFkX4c/Y6sYTTc5kyZYo7c1+bDkC1PSKtt87Ke1NLine truncated
"text/plain": [
"<Figure size 1100x850 with 1 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"from matplotlib.patches import Ellipse\n",
"\n",
"# the tight Kiddieland pile = rides loading hard-negative on PC2\n",
"kiddie_core = loadings.index[loadings[\"PC2\"] < -0.13]\n",
"kx, ky = loadings.loc[kiddie_core, \"PC2\"], loadings.loc[kiddie_core, \"PC3\"]\n",
"\n",
"fig, ax = plt.subplots(figsize=(11, 8.5))\n",
"cat_colors = {\"Family\": \"tab:blue\", \"Kiddie\": \"tab:orange\", \"Thrill\": \"tab:red\"}\n",
"for c, g in annotated.dropna(subset=[\"category\"]).groupby(\"category\"):\n",
" ax.scatter(g[\"PC2\"], g[\"PC3\"], label=c, color=cat_colors.get(c, \"grey\"),\n",
" s=80, alpha=0.75, edgecolor=\"white\", zorder=3)\n",
"\n",
"# one callout for the dense kiddie pile, instead of 13 overlapping labels\n",
"w, h = (kx.max() - kx.min()) + 0.07, (ky.max() - ky.min()) + 0.07\n",
"ax.add_patch(Ellipse((kx.mean(), ky.mean()), width=w, height=h, facecolor=\"tab:orange\",\n",
" alpha=0.08, edgecolor=\"tab:orange\", ls=\"--\", lw=1.5, zorder=1))\n",
"ax.annotate(\"Kiddieland\\n(tight cluster of kiddie rides)\", (kx.mean(), ky.mean() - h / 2),\n",
" ha=\"center\", va=\"top\", fontsize=9, color=\"tab:orange\", fontweight=\"bold\")\n",
"\n",
"# label only the rides that carry the story\n",
"plain = {\"giant_wheel\": (6, 5), \"grand_carousel\": (8, -3), \"impulse\": (7, -3),\n",
" \"bumper_cars\": (7, 4), \"motor_boats\": (7, 0), \"paratrooper\": (7, 0),\n",
" \"tea_cups\": (7, 0), \"tilt-a-whirl\": (7, -2), \"roto_jets\": (-42, 6)}\n",
"for ride, (dx, dy) in plain.items():\n",
" p = annotated.loc[ride]\n",
" ax.annotate(ride, (p[\"PC2\"], p[\"PC3\"]), fontsize=8.5,\n",
" xytext=(dx, dy), textcoords=\"offset points\", zorder=4)\n",
"\n",
"# the two Family-tagged rides that behave like kiddie rides \u2014 pull labels out with arrows\n",
"misfiled = {\"s&g_carousel\": (-70, -30), \"kiddie_himalaya\": (-78, 20)}\n",
"for ride, (dx, dy) in misfiled.items():\n",
" p = annotated.loc[ride]\n",
" ax.annotate(ride + \" (tagged Family)\", (p[\"PC2\"], p[\"PC3\"]), fontsize=8.5,\n",
" color=\"tab:blue\", xytext=(dx, dy), textcoords=\"offset points\",\n",
" arrowprops=dict(arrowstyle=\"->\", lw=0.7, color=\"tab:blue\"), zorder=5)\n",
"\n",
"ax.axhline(0, color=\"grey\", lw=0.5)\n",
"ax.axvline(0, color=\"grey\", lw=0.5)\n",
"ax.set_xlabel(\"PC2 Kiddieland corner \u2190 \u2192 main midway\")\n",
"ax.set_ylabel(\"PC3 evening crowd \u2190 \u2192 daytime / school-group daypart\")\n",
"ax.set_title(\"Knoebels rides in PC2 \u00d7 PC3 space\\n\"\n",
" \"(position comes from queue patterns alone; color = the park's category tag)\")\n",
"ax.legend(title=\"Category tag\", loc=\"lower left\")\n",
"plt.tight_layout()"
]
},
{
"cell_type": "markdown",
"id": "73b2a3f7",
"metadata": {},
"source": [
"## What I took away\n",
"\n",
"- **PCA found real structure with zero domain input.** Three axes \u2014 busyness,\n",
" geography, time of day \u2014 fell out of nothing but wait times, and each one maps\n",
" to something you could verify against the park independently (the average wait,\n",
" the GPS map, the operating schedule).\n",
"- **Standardizing matters twice.** Once before PCA (so small rides aren't\n",
" drowned out), and again before k-means on coordinates (so the clustering is\n",
" reproducible).\n",
"- **The model's labels beat the official ones.** Where the category tag and the\n",
" queue behavior disagreed (`giant_wheel`, `s&g_carousel`), the behavior was the\n",
" truer description of how a ride is actually *used*.\n",
"- **Be willing to kill a pretty hypothesis.** The PC3 \"throughput\" story was\n",
" tidy and wrong; the data pointed at time-of-day instead.\n",
"\n",
"Skipping ahead in the textbook to chase a dataset I cared about turned out to be\n",
"a great way to learn this."
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"name": "python"
}
},
"nbformat": 4,
"nbformat_minor": 5
}