{"metadata":{"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"},"language_info":{"name":"python","version":"3.12.13","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"}},"nbformat_minor":4,"nbformat":4,"cells":[{"id":"6e1286b0-5772-4c66-88dc-018955b156ff","cell_type":"markdown","source":"# Global Earthquake Tsunami Risk Assessment\n### Kaggle Dataset | Tsunami Classification\n\n**Goal:** predict whether a recorded earthquake triggers a tsunami, using seismic and location features from the USGS-style event catalog.\n\n**Target:** `tsunami` (0 or 1). About 6 percent of events are positive, so this is an imbalanced classification problem.\n\n**Evaluation metric:** ROC-AUC as the primary metric, with precision, recall, and F1 tracked alongside it, since accuracy alone is misleading on an imbalanced target.\n\n**Approach:**\n1. Setup\n2. Data loading and cleaning\n3. Exploratory data analysis\n4. Feature engineering\n5. Dimensionality reduction (PCA)\n6. Clustering (K-Means)\n7. Baseline classification models (XGBoost, LightGBM)\n8. Hyperparameter tuning with Optuna\n9. Final evaluation","metadata":{}},{"id":"ed2d918c-1160-4898-82d4-d74825250b32","cell_type":"markdown","source":"---\n## 1. Setup","metadata":{}},{"id":"184667e4-baa2-45e6-9068-2261e6bfd37f","cell_type":"code","source":"import warnings\nimport numpy as np\nimport pandas as pd\nimport matplotlib.pyplot as plt\nimport seaborn as sns\n\nfrom sklearn.model_selection import train_test_split\nfrom sklearn.preprocessing import LabelEncoder, StandardScaler\nfrom sklearn.decomposition import PCA\nfrom sklearn.cluster import KMeans\nfrom sklearn.metrics import (\n    accuracy_score, precision_score, recall_score, f1_score,\n    roc_auc_score, roc_curve, confusion_matrix, classification_report\n)\n\nimport xgboost as xgb\nimport lightgbm as lgb\nimport optuna\n\nwarnings.filterwarnings(\"ignore\")\noptuna.logging.set_verbosity(optuna.logging.WARNING)\n\nSEED = 42\nsns.set_style(\"whitegrid\")\nplt.rcParams[\"figure.figsize\"] = (10, 6)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-08-26T11:29:20.170616Z","iopub.execute_input":"2026-08-26T11:29:20.170938Z","iopub.status.idle":"2026-08-26T11:29:20.177780Z","shell.execute_reply.started":"2026-08-26T11:29:20.170912Z","shell.execute_reply":"2026-08-26T11:29:20.176740Z"}},"outputs":[],"execution_count":null},{"id":"f2a1cf2c-5025-4a05-aa44-ad71c129c1af","cell_type":"markdown","source":"---\n## 2. Data Loading and Cleaning","metadata":{}},{"id":"db2f2c70-fbeb-4e47-bcae-7d9f460fcb87","cell_type":"markdown","source":"### 2.1 Load and inspect","metadata":{}},{"id":"4bf742d9-61e2-446b-9f92-12e279017e2b","cell_type":"code","source":"df = pd.read_csv(\"/kaggle/input/datasets/shreyasur965/recent-earthquakes/earthquakes.csv\")\nprint(df.shape)\ndf.head()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-08-26T11:29:20.179256Z","iopub.execute_input":"2026-08-26T11:29:20.179586Z","iopub.status.idle":"2026-08-26T11:29:20.249668Z","shell.execute_reply.started":"2026-08-26T11:29:20.179535Z","shell.execute_reply":"2026-08-26T11:29:20.248835Z"}},"outputs":[],"execution_count":null},{"id":"332aeb7f-4d92-44b0-86bf-f5cc215c35da","cell_type":"code","source":"df.info()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-08-26T11:29:20.251253Z","iopub.execute_input":"2026-08-26T11:29:20.251654Z","iopub.status.idle":"2026-08-26T11:29:20.271306Z","shell.execute_reply.started":"2026-08-26T11:29:20.251620Z","shell.execute_reply":"2026-08-26T11:29:20.270319Z"}},"outputs":[],"execution_count":null},{"id":"c4db14ff-9ce1-4d63-9136-48f08e7fe14c","cell_type":"code","source":"missing = df.isnull().sum()\nmissing[missing > 0].sort_values(ascending=False)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-08-26T11:29:20.272361Z","iopub.execute_input":"2026-08-26T11:29:20.273288Z","iopub.status.idle":"2026-08-26T11:29:20.294055Z","shell.execute_reply.started":"2026-08-26T11:29:20.273243Z","shell.execute_reply":"2026-08-26T11:29:20.293122Z"}},"outputs":[],"execution_count":null},{"id":"f3811c8b-b258-414d-a341-4da1da3169c5","cell_type":"markdown","source":"alert, continent, country, subnational, city, postcode all have gaps, everything else is complete\n\n=> alert is missing when USGS did not issue a PAGER alert, so fill with none instead of dropping rows\n=> continent/country gaps get filled as Unknown so no rows are lost","metadata":{}},{"id":"77d7bfd2-ab16-4ac5-b823-4a5afca03a72","cell_type":"markdown","source":"### 2.2 Drop low-signal columns and fill gaps","metadata":{}},{"id":"d7eb74b2-8349-4595-892a-b94b697589ce","cell_type":"code","source":"DROP_COLS = [\n    \"url\", \"detailUrl\", \"title\", \"code\", \"ids\", \"types\", \"what3words\",\n    \"locationDetails\", \"postcode\", \"placeOnly\", \"location\", \"subnational\",\n    \"city\", \"locality\", \"net\", \"time\", \"updated\", \"timezone\"\n]\n\ncleaned_df = df.drop(columns=[c for c in DROP_COLS if c in df.columns]).copy()\n\ncleaned_df[\"alert\"] = cleaned_df[\"alert\"].fillna(\"none\")\ncleaned_df[\"continent\"] = cleaned_df[\"continent\"].fillna(\"Unknown\")\ncleaned_df[\"country\"] = cleaned_df[\"country\"].fillna(\"Unknown\")\n\ncleaned_df[\"date\"] = pd.to_datetime(cleaned_df[\"date\"])\ncleaned_df[\"year\"] = cleaned_df[\"date\"].dt.year\ncleaned_df[\"month\"] = cleaned_df[\"date\"].dt.month\ncleaned_df[\"day_of_week\"] = cleaned_df[\"date\"].dt.dayofweek\ncleaned_df[\"hour\"] = cleaned_df[\"date\"].dt.hour\n\ncleaned_df = cleaned_df.dropna(subset=[\"magnitude\", \"depth\", \"latitude\", \"longitude\"])\nprint(\"Cleaned shape:\", cleaned_df.shape)\ncleaned_df.isnull().sum().sum()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-08-26T11:29:20.295235Z","iopub.execute_input":"2026-08-26T11:29:20.295578Z","iopub.status.idle":"2026-08-26T11:29:20.330418Z","shell.execute_reply.started":"2026-08-26T11:29:20.295528Z","shell.execute_reply":"2026-08-26T11:29:20.329597Z"}},"outputs":[],"execution_count":null},{"id":"5b2f467c-6103-48b1-ae50-8c1550b1c523","cell_type":"markdown","source":"dropped columns are IDs, URLs, or free text with little predictive value\n\nzero nulls left in cleaned_df","metadata":{}},{"id":"7247fef2-01e9-484f-9f3e-be320f684140","cell_type":"markdown","source":"---\n## 3. Exploratory Data Analysis","metadata":{}},{"id":"96ede42b-7ff2-45db-97c1-1f9b311c37ba","cell_type":"markdown","source":"### 3.1 Magnitude distribution","metadata":{}},{"id":"d0798b93-fedc-4ae3-9d7b-9c69d757d671","cell_type":"code","source":"plt.figure(figsize=(10, 6))\nsns.histplot(cleaned_df[\"magnitude\"], bins=30, kde=True, color=\"steelblue\")\nplt.title(\"Distribution of Earthquake Magnitudes\")\nplt.xlabel(\"Magnitude\")\nplt.ylabel(\"Frequency\")\nplt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-08-26T11:29:20.332813Z","iopub.execute_input":"2026-08-26T11:29:20.333086Z","iopub.status.idle":"2026-08-26T11:29:20.675551Z","shell.execute_reply.started":"2026-08-26T11:29:20.333060Z","shell.execute_reply":"2026-08-26T11:29:20.674750Z"}},"outputs":[],"execution_count":null},{"id":"cd78f464-231c-4ed3-a6b0-18d3b77b2c6c","cell_type":"markdown","source":"most events cluster between magnitude 4 and 5, with a long right tail toward the rare, larger events","metadata":{}},{"id":"545f6b16-e769-4214-840b-1582590b83a7","cell_type":"markdown","source":"### 3.2 Depth distribution","metadata":{}},{"id":"540cfbd9-33b9-4e5b-ba7f-65e14ad5b124","cell_type":"code","source":"plt.figure(figsize=(10, 6))\nsns.histplot(cleaned_df[\"depth\"], bins=30, kde=True, color=\"seagreen\")\nplt.title(\"Distribution of Earthquake Depths\")\nplt.xlabel(\"Depth (km)\")\nplt.ylabel(\"Frequency\")\nplt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-08-26T11:29:20.676707Z","iopub.execute_input":"2026-08-26T11:29:20.677038Z","iopub.status.idle":"2026-08-26T11:29:20.910379Z","shell.execute_reply.started":"2026-08-26T11:29:20.677005Z","shell.execute_reply":"2026-08-26T11:29:20.909661Z"}},"outputs":[],"execution_count":null},{"id":"6490a67b-d436-4d0c-9ed3-4ab2b751e386","cell_type":"markdown","source":"depth is also right-skewed, most events are shallow (under 50 km) with a long tail of deep events","metadata":{}},{"id":"b547819e-b82d-4161-a1ca-1c681f7c87e4","cell_type":"markdown","source":"### 3.3 Depth vs magnitude, split by tsunami","metadata":{}},{"id":"2ea450a5-0ac2-4a4b-b5ef-c2be9043b59d","cell_type":"code","source":"plt.figure(figsize=(10, 6))\nsns.scatterplot(x=\"depth\", y=\"magnitude\", hue=\"tsunami\", data=cleaned_df,\n                 palette={0: \"lightgray\", 1: \"crimson\"}, alpha=0.7)\nplt.title(\"Depth vs Magnitude, Colored by Tsunami Occurrence\")\nplt.xlabel(\"Depth (km)\")\nplt.ylabel(\"Magnitude\")\nplt.legend(title=\"Tsunami\")\nplt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-08-26T11:29:20.911402Z","iopub.execute_input":"2026-08-26T11:29:20.912375Z","iopub.status.idle":"2026-08-26T11:29:21.152790Z","shell.execute_reply.started":"2026-08-26T11:29:20.912314Z","shell.execute_reply":"2026-08-26T11:29:21.151913Z"}},"outputs":[],"execution_count":null},{"id":"0bb69948-14f3-4050-88bd-4db45a513cb3","cell_type":"markdown","source":"tsunami-triggering events cluster at higher magnitude and mostly shallow to moderate depth\n\n=> magnitude and depth are likely strong predictors","metadata":{}},{"id":"3bd85da5-9b2e-4604-be0e-bcc35a06e8d0","cell_type":"markdown","source":"### 3.4 Feature correlation with the tsunami target\n\nThe heatmap in the next section shows how features relate to each other, this plot isolates just the correlation of each numeric feature with tsunami itself and sorts it, so the strongest single predictors are easy to read off directly.","metadata":{}},{"id":"5e4d7a07-cf65-49a4-ab7b-dae4b5084f3b","cell_type":"code","source":"NUMERIC_COLS = [\"magnitude\", \"felt\", \"cdi\", \"mmi\", \"sig\", \"nst\", \"dmin\", \"rms\",\n                \"gap\", \"depth\", \"latitude\", \"longitude\", \"distanceKM\"]\n\ntarget_corr = cleaned_df[NUMERIC_COLS + [\"tsunami\"]].corr()[\"tsunami\"].drop(\"tsunami\")\ntarget_corr = target_corr.sort_values()\n\nplt.figure(figsize=(9, 7))\ncolors = [\"crimson\" if v > 0 else \"steelblue\" for v in target_corr.values]\nsns.barplot(x=target_corr.values, y=target_corr.index, palette=colors)\nplt.title(\"Correlation of Each Feature with Tsunami Occurrence\")\nplt.xlabel(\"Correlation with tsunami\")\nplt.axvline(0, color=\"black\", lw=0.8)\nplt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-08-26T11:29:21.153952Z","iopub.execute_input":"2026-08-26T11:29:21.154291Z","iopub.status.idle":"2026-08-26T11:29:21.440264Z","shell.execute_reply.started":"2026-08-26T11:29:21.154257Z","shell.execute_reply":"2026-08-26T11:29:21.439450Z"}},"outputs":[],"execution_count":null},{"id":"95c78785-6344-463e-9f7f-12eb949cb7ca","cell_type":"markdown","source":"sig and magnitude carry the strongest positive correlation with tsunami, gap and nst are close to zero and add little on their own","metadata":{}},{"id":"114304ae-b1e4-4aea-bdbe-6ab41dd1cf91","cell_type":"markdown","source":"### 3.5 Correlation heatmap","metadata":{}},{"id":"d7d2f061-8f08-426a-8449-75f0de78e87b","cell_type":"code","source":"plt.figure(figsize=(13, 10))\nsns.heatmap(cleaned_df[NUMERIC_COLS + [\"tsunami\"]].corr(), annot=True, fmt=\".2f\", cmap=\"coolwarm\")\nplt.title(\"Correlation Heatmap of Numerical Features\")\nplt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-08-26T11:29:21.441407Z","iopub.execute_input":"2026-08-26T11:29:21.442245Z","iopub.status.idle":"2026-08-26T11:29:22.008840Z","shell.execute_reply.started":"2026-08-26T11:29:21.442205Z","shell.execute_reply":"2026-08-26T11:29:22.007700Z"}},"outputs":[],"execution_count":null},{"id":"c7314d59-088e-4983-9479-ca287c6ae991","cell_type":"markdown","source":"sig, cdi, mmi, and felt all move together since sig is partly derived from them, worth keeping in mind so the model is not just relearning one signal five times over","metadata":{}},{"id":"b9baa888-1745-4dc4-b5e8-be02c3591ab9","cell_type":"markdown","source":"### 3.6 Pairwise relationships between key seismic measurements\n\nA pairplot puts every pair of the strongest candidate features against each other in one grid, faster than building each scatter plot one at a time.","metadata":{}},{"id":"27ed6781-90fd-4e6a-83e2-056d98b02a9c","cell_type":"code","source":"pair_cols = [\"magnitude\", \"depth\", \"sig\", \"cdi\", \"mmi\", \"tsunami\"]\nsns.pairplot(cleaned_df[pair_cols], hue=\"tsunami\", palette={0: \"lightgray\", 1: \"crimson\"},\n             diag_kind=\"kde\", plot_kws={\"alpha\": 0.6, \"s\": 20})\nplt.suptitle(\"Pairplot of Key Seismic Measurements by Tsunami Occurrence\", y=1.02)\nplt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-08-26T11:29:22.010426Z","iopub.execute_input":"2026-08-26T11:29:22.011179Z","iopub.status.idle":"2026-08-26T11:29:27.023297Z","shell.execute_reply.started":"2026-08-26T11:29:22.011119Z","shell.execute_reply":"2026-08-26T11:29:27.022449Z"}},"outputs":[],"execution_count":null},{"id":"74b7e9a2-7650-4ac4-af2c-ede382d38513","cell_type":"markdown","source":"tsunami points separate out most clearly on the magnitude and sig panels, and barely separate at all on cdi/mmi alone","metadata":{}},{"id":"53e9b4e4-6bae-4dae-aa4b-2a9f037cad0c","cell_type":"markdown","source":"### 3.7 Target balance","metadata":{}},{"id":"d140d515-f10f-42cd-bb52-a18e0c14222d","cell_type":"code","source":"tsunami_counts = cleaned_df[\"tsunami\"].value_counts()\nprint(tsunami_counts)\nprint(\"Positive class share: {:.2f} percent\".format(100 * tsunami_counts[1] / tsunami_counts.sum()))\n\nplt.figure(figsize=(7, 5))\nsns.barplot(x=[\"No Tsunami\", \"Tsunami\"], y=tsunami_counts.values, palette=[\"lightblue\", \"salmon\"])\nplt.title(\"Tsunami Occurrence Counts\")\nplt.ylabel(\"Number of Events\")\nplt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-08-26T11:29:27.024529Z","iopub.execute_input":"2026-08-26T11:29:27.024866Z","iopub.status.idle":"2026-08-26T11:29:27.163307Z","shell.execute_reply.started":"2026-08-26T11:29:27.024840Z","shell.execute_reply":"2026-08-26T11:29:27.162414Z"}},"outputs":[],"execution_count":null},{"id":"e4e28458-ad5e-48c7-a5c3-39ba6cf2a810","cell_type":"markdown","source":"about 6 percent positive class\n\n=> use a stratified split and class weighting later, plain accuracy would be misleading here","metadata":{}},{"id":"3949ed45-0333-4551-9ed7-ef52e31b6719","cell_type":"markdown","source":"### 3.8 Magnitude and significance by tsunami occurrence","metadata":{}},{"id":"e0622685-9e4a-4521-9b5a-00ae3423c911","cell_type":"code","source":"fig, axes = plt.subplots(1, 2, figsize=(13, 6))\nsns.boxplot(x=\"tsunami\", y=\"magnitude\", data=cleaned_df, ax=axes[0], palette=[\"lightblue\", \"salmon\"])\naxes[0].set_title(\"Magnitude by Tsunami Occurrence\")\naxes[0].set_xlabel(\"Tsunami\")\n\nsns.boxplot(x=\"tsunami\", y=\"sig\", data=cleaned_df, ax=axes[1], palette=[\"lightblue\", \"salmon\"])\naxes[1].set_title(\"Significance Score by Tsunami Occurrence\")\naxes[1].set_xlabel(\"Tsunami\")\n\nplt.tight_layout()\nplt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-08-26T11:29:27.164386Z","iopub.execute_input":"2026-08-26T11:29:27.164756Z","iopub.status.idle":"2026-08-26T11:29:27.469467Z","shell.execute_reply.started":"2026-08-26T11:29:27.164723Z","shell.execute_reply":"2026-08-26T11:29:27.468726Z"}},"outputs":[],"execution_count":null},{"id":"a5459d38-25ce-4c93-96d0-a01fdbd6ad3e","cell_type":"markdown","source":"both distributions shift up clearly for tsunami events, with sig showing an even sharper split than magnitude alone","metadata":{}},{"id":"d25e5134-5dcb-4b94-9901-855bd59eceac","cell_type":"markdown","source":"### 3.9 Magnitude by continent","metadata":{}},{"id":"02bd1c19-6197-4270-83a6-cb8e44b3b50c","cell_type":"code","source":"plt.figure(figsize=(11, 7))\nsns.boxplot(x=\"continent\", y=\"magnitude\", data=cleaned_df)\nplt.title(\"Earthquake Magnitude by Continent\")\nplt.xlabel(\"Continent\")\nplt.ylabel(\"Magnitude\")\nplt.xticks(rotation=40, ha=\"right\")\nplt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-08-26T11:29:27.470514Z","iopub.execute_input":"2026-08-26T11:29:27.470811Z","iopub.status.idle":"2026-08-26T11:29:27.679240Z","shell.execute_reply.started":"2026-08-26T11:29:27.470788Z","shell.execute_reply":"2026-08-26T11:29:27.678366Z"}},"outputs":[],"execution_count":null},{"id":"e46d5499-e585-4ae0-aa98-ee52d758709e","cell_type":"markdown","source":"### 3.10 Geographic distribution","metadata":{}},{"id":"ce9b066c-4961-4107-a6d8-f76e0b04b0f9","cell_type":"code","source":"plt.figure(figsize=(11, 7))\nsns.scatterplot(x=\"longitude\", y=\"latitude\", hue=\"tsunami\", size=\"magnitude\",\n                 data=cleaned_df, palette={0: \"lightgray\", 1: \"crimson\"}, alpha=0.75)\nplt.title(\"Geographic Distribution of Earthquakes and Tsunami Occurrence\")\nplt.xlabel(\"Longitude\")\nplt.ylabel(\"Latitude\")\nplt.legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.)\nplt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-08-26T11:29:27.683069Z","iopub.execute_input":"2026-08-26T11:29:27.683377Z","iopub.status.idle":"2026-08-26T11:29:28.156419Z","shell.execute_reply.started":"2026-08-26T11:29:27.683351Z","shell.execute_reply":"2026-08-26T11:29:28.155626Z"}},"outputs":[],"execution_count":null},{"id":"18e38dbd-e43d-48cf-a0ed-0fc9d8b4a874","cell_type":"markdown","source":"tsunami events sit mostly along the Pacific coastlines, matching the known Ring of Fire subduction zones\n\n=> latitude/longitude carry real signal even without an explicit fault-line feature","metadata":{}},{"id":"74b7b2d8-d84e-4caf-b462-1139a593d2e5","cell_type":"markdown","source":"### 3.11 Alert level vs tsunami","metadata":{}},{"id":"cf49f0cf-b849-4610-9124-bd6f974e1245","cell_type":"code","source":"tsunami_alert = cleaned_df.groupby([\"alert\", \"tsunami\"]).size().unstack(fill_value=0)\ntsunami_alert.plot(kind=\"bar\", stacked=True, figsize=(9, 6), color=[\"lightblue\", \"salmon\"])\nplt.title(\"Tsunami Occurrence by Alert Level\")\nplt.xlabel(\"Alert Level\")\nplt.ylabel(\"Number of Events\")\nplt.legend([\"No Tsunami\", \"Tsunami\"])\nplt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-08-26T11:29:28.157519Z","iopub.execute_input":"2026-08-26T11:29:28.157874Z","iopub.status.idle":"2026-08-26T11:29:28.410317Z","shell.execute_reply.started":"2026-08-26T11:29:28.157841Z","shell.execute_reply":"2026-08-26T11:29:28.409596Z"}},"outputs":[],"execution_count":null},{"id":"cc51f9c9-cc1f-4e83-b59c-a532ec302158","cell_type":"markdown","source":"orange and red alerts have a much higher tsunami rate than green, as expected since alert level is itself partly driven by impact severity","metadata":{}},{"id":"3ab8e943-d77b-48a7-bfa3-e4b617beeeed","cell_type":"markdown","source":"### 3.12 Earthquakes over time","metadata":{}},{"id":"14bb1989-009d-40d2-8bc0-d32bb15b1abb","cell_type":"code","source":"monthly_counts = cleaned_df.set_index(\"date\").resample(\"ME\").size()\nmonthly_tsunami = cleaned_df[cleaned_df[\"tsunami\"] == 1].set_index(\"date\").resample(\"ME\").size()\n\nfig, ax = plt.subplots(figsize=(14, 6))\nmonthly_counts.plot(ax=ax, color=\"steelblue\", label=\"All Earthquakes\")\nmonthly_tsunami.plot(ax=ax, color=\"crimson\", label=\"Tsunami Events\")\nax.set_title(\"Earthquake and Tsunami Counts Over Time\")\nax.set_xlabel(\"Date\")\nax.set_ylabel(\"Number of Events per Month\")\nax.legend()\nplt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-08-26T11:29:28.411375Z","iopub.execute_input":"2026-08-26T11:29:28.411686Z","iopub.status.idle":"2026-08-26T11:29:28.646238Z","shell.execute_reply.started":"2026-08-26T11:29:28.411652Z","shell.execute_reply":"2026-08-26T11:29:28.645338Z"}},"outputs":[],"execution_count":null},{"id":"dffffab8-cb24-4f07-93f5-a9ea222b325c","cell_type":"markdown","source":"tsunami counts track total earthquake activity rather than showing an independent seasonal pattern of their own","metadata":{}},{"id":"4536b938-9739-4c2b-a2e0-7beffa299da6","cell_type":"markdown","source":"---\n## 4. Feature Engineering\n\nCategorical columns are label-encoded and combined with the numeric columns into one feature matrix, then scaled. Scaling matters here because PCA and K-Means both use distances between points, and unscaled columns like sig (hundreds) would dominate over columns like magnitude (single digits).\n\n| Group | Columns |\n|---|---|\n| Numeric | magnitude, felt, cdi, mmi, sig, nst, dmin, rms, gap, depth, latitude, longitude, distanceKM |\n| Calendar | year, month, day_of_week, hour |\n| Categorical (encoded) | alert, magType, continent, status, type |","metadata":{}},{"id":"3796e7b9-1095-4732-bd04-d06a477c7acb","cell_type":"code","source":"FEATURE_COLS = [\"magnitude\", \"felt\", \"cdi\", \"mmi\", \"sig\", \"nst\", \"dmin\", \"rms\", \"gap\",\n                \"depth\", \"latitude\", \"longitude\", \"distanceKM\", \"year\", \"month\",\n                \"day_of_week\", \"hour\"]\n\nCAT_COLS = [\"alert\", \"magType\", \"continent\", \"status\", \"type\"]\n\nmodel_df = cleaned_df[FEATURE_COLS + CAT_COLS + [\"tsunami\"]].copy()\n\nencoders = {}\nfor c in CAT_COLS:\n    le = LabelEncoder()\n    model_df[c] = le.fit_transform(model_df[c].astype(str))\n    encoders[c] = le\n\nX = model_df.drop(columns=[\"tsunami\"])\ny = model_df[\"tsunami\"]\n\nscaler = StandardScaler()\nX_scaled = scaler.fit_transform(X)\nprint(X.shape, y.shape)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-08-26T11:29:28.647341Z","iopub.execute_input":"2026-08-26T11:29:28.647688Z","iopub.status.idle":"2026-08-26T11:29:28.666083Z","shell.execute_reply.started":"2026-08-26T11:29:28.647662Z","shell.execute_reply":"2026-08-26T11:29:28.665150Z"}},"outputs":[],"execution_count":null},{"id":"7b7dd63c-2bb5-433b-a05a-6e5a725f4a2a","cell_type":"markdown","source":"---\n## 5. Dimensionality Reduction\n\n22 features cannot be plotted directly, so PCA is used to project them down to 2 dimensions and check by eye whether tsunami and non-tsunami events form separable groups.","metadata":{}},{"id":"07caf688-1635-4982-bdc0-d073c706a7d4","cell_type":"markdown","source":"(might remove it actually)","metadata":{}},{"id":"7d906cd8-2e0f-4f86-a0f4-5c4b1e6dab03","cell_type":"markdown","source":"### 5.1 PCA \n\nPCA finds the directions (principal components) that capture the most variance in the data and projects onto the top two. It is linear and fast, and the axes have a clear meaning: each is a weighted combination of the original features.","metadata":{}},{"id":"238c604f-4eeb-489f-9b02-27a47921a94d","cell_type":"code","source":"pca = PCA(n_components=2, random_state=SEED)\nX_pca = pca.fit_transform(X_scaled)\nprint(\"Explained variance ratio:\", pca.explained_variance_ratio_)\n\nplt.figure(figsize=(9, 7))\nsns.scatterplot(x=X_pca[:, 0], y=X_pca[:, 1], hue=y, palette={0: \"lightgray\", 1: \"crimson\"}, alpha=0.7)\nplt.title(\"PCA Projection of Seismic Features\")\nplt.xlabel(\"PC1\")\nplt.ylabel(\"PC2\")\nplt.legend(title=\"Tsunami\")\nplt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-08-26T11:38:26.434682Z","iopub.execute_input":"2026-08-26T11:38:26.435024Z","iopub.status.idle":"2026-08-26T11:38:26.655150Z","shell.execute_reply.started":"2026-08-26T11:38:26.434997Z","shell.execute_reply":"2026-08-26T11:38:26.654317Z"}},"outputs":[],"execution_count":null},{"id":"12913efe-8b68-4717-aab8-49a373360ab2","cell_type":"markdown","source":"the two components only capture a modest share of total variance, so PCA gives partial but not clean separation","metadata":{}},{"id":"0db41578-06c1-4538-9647-530fbd7ce63b","cell_type":"markdown","source":"---\n## 6. Clustering\n\nClustering is unsupervised: the tsunami label is not used. The question here is whether earthquakes naturally group into risk profiles based on their seismic features alone.","metadata":{}},{"id":"e42adb18-98d8-492f-96f2-4d04f77638e1","cell_type":"markdown","source":"### 6.1 K-Means\n\nK-Means requires the number of clusters (k) to be chosen up front. The elbow plot below runs K-Means for a range of k values and plots inertia (within-cluster distance), so the leveling-off point suggests a reasonable k.","metadata":{}},{"id":"11566a8e-42ea-4f3a-99f8-2f985f4cc265","cell_type":"code","source":"inertias = []\nk_range = range(2, 11)\nfor k in k_range:\n    km = KMeans(n_clusters=k, random_state=SEED, n_init=10)\n    km.fit(X_scaled)\n    inertias.append(km.inertia_)\n\nplt.figure(figsize=(9, 6))\nplt.plot(list(k_range), inertias, marker=\"o\")\nplt.title(\"K-Means Elbow Plot\")\nplt.xlabel(\"Number of Clusters\")\nplt.ylabel(\"Inertia\")\nplt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-08-26T11:29:28.901505Z","iopub.execute_input":"2026-08-26T11:29:28.901828Z","iopub.status.idle":"2026-08-26T11:29:29.372193Z","shell.execute_reply.started":"2026-08-26T11:29:28.901805Z","shell.execute_reply":"2026-08-26T11:29:29.371312Z"}},"outputs":[],"execution_count":null},{"id":"c98546f6-49e9-4b32-b30c-33ba9b38b777","cell_type":"code","source":"kmeans = KMeans(n_clusters=4, random_state=SEED, n_init=10)\ncluster_labels = kmeans.fit_predict(X_scaled)\n\nplt.figure(figsize=(9, 7))\nsns.scatterplot(x=X_pca[:, 0], y=X_pca[:, 1], hue=cluster_labels, palette=\"tab10\", alpha=0.7)\nplt.title(\"K-Means Clusters Projected on PCA Components\")\nplt.xlabel(\"PC1\")\nplt.ylabel(\"PC2\")\nplt.legend(title=\"Cluster\")\nplt.show()\n\npd.crosstab(cluster_labels, y, rownames=[\"Cluster\"], colnames=[\"Tsunami\"])","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-08-26T11:29:29.373351Z","iopub.execute_input":"2026-08-26T11:29:29.373685Z","iopub.status.idle":"2026-08-26T11:29:29.638203Z","shell.execute_reply.started":"2026-08-26T11:29:29.373662Z","shell.execute_reply":"2026-08-26T11:29:29.637507Z"}},"outputs":[],"execution_count":null},{"id":"a4f9496d-8c82-4dfe-b821-fe7fafe33eea","cell_type":"markdown","source":"one or two clusters carry a disproportionate share of tsunami events, so K-Means partially recovers the risk structure without ever seeing the label","metadata":{}},{"id":"6699e868-0f71-4935-b73b-06ae6d3078a1","cell_type":"markdown","source":"---\n## 7. Baseline Classification Models\n\nThe target is imbalanced (about 6 percent positive), so the split is stratified to keep that ratio in both train and test, and each model is given a way to weight the minority class more heavily during training.","metadata":{}},{"id":"726454ac-7942-477c-b063-e358a4263a7f","cell_type":"code","source":"X_train, X_test, y_train, y_test = train_test_split(\n    X, y, test_size=0.2, random_state=SEED, stratify=y\n)\nprint(\"Train shape:\", X_train.shape, \" Test shape:\", X_test.shape)\nprint(\"Train positive rate:\", round(y_train.mean(), 4))\nprint(\"Test positive rate:\", round(y_test.mean(), 4))","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-08-26T11:29:29.639166Z","iopub.execute_input":"2026-08-26T11:29:29.639412Z","iopub.status.idle":"2026-08-26T11:29:29.650069Z","shell.execute_reply.started":"2026-08-26T11:29:29.639390Z","shell.execute_reply":"2026-08-26T11:29:29.649088Z"}},"outputs":[],"execution_count":null},{"id":"5573b2fe-ab1d-46c2-bf44-0c94464327d6","cell_type":"markdown","source":"XGBoost and LightGBM are compared as baseline gradient boosting models, both run with comparable settings so the results are a fair comparison, not a tuned one, tuning comes later.","metadata":{}},{"id":"6d60bd01-9750-4fcd-bf6f-be3af2f7b305","cell_type":"code","source":"def evaluate_model(name, model, X_test, y_test):\n    preds = model.predict(X_test)\n    proba = model.predict_proba(X_test)[:, 1]\n    return {\n        \"model\": name,\n        \"accuracy\": accuracy_score(y_test, preds),\n        \"precision\": precision_score(y_test, preds, zero_division=0),\n        \"recall\": recall_score(y_test, preds, zero_division=0),\n        \"f1\": f1_score(y_test, preds, zero_division=0),\n        \"roc_auc\": roc_auc_score(y_test, proba)\n    }, proba\n\nresults = []\nproba_by_model = {}\nscale_pos_weight = (y_train == 0).sum() / (y_train == 1).sum()\n\nxgb_model = xgb.XGBClassifier(\n    n_estimators=300, max_depth=5, learning_rate=0.05,\n    scale_pos_weight=scale_pos_weight, eval_metric=\"logloss\", random_state=SEED\n)\nxgb_model.fit(X_train, y_train)\nm, p = evaluate_model(\"XGBoost\", xgb_model, X_test, y_test)\nresults.append(m); proba_by_model[\"XGBoost\"] = p\n\nlgb_model = lgb.LGBMClassifier(\n    n_estimators=300, max_depth=5, learning_rate=0.05,\n    class_weight=\"balanced\", random_state=SEED, verbose=-1\n)\nlgb_model.fit(X_train, y_train)\nm, p = evaluate_model(\"LightGBM\", lgb_model, X_test, y_test)\nresults.append(m); proba_by_model[\"LightGBM\"] = p\n\nresults_df = pd.DataFrame(results).set_index(\"model\").round(4)\nresults_df","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-08-26T11:29:29.651125Z","iopub.execute_input":"2026-08-26T11:29:29.651755Z","iopub.status.idle":"2026-08-26T11:29:30.180924Z","shell.execute_reply.started":"2026-08-26T11:29:29.651729Z","shell.execute_reply":"2026-08-26T11:29:30.180080Z"}},"outputs":[],"execution_count":null},{"id":"34cabe1d-f347-46bd-b0fd-da4af9fdcc43","cell_type":"code","source":"plt.figure(figsize=(9, 7))\nfor name, proba in proba_by_model.items():\n    fpr, tpr, _ = roc_curve(y_test, proba)\n    auc = roc_auc_score(y_test, proba)\n    plt.plot(fpr, tpr, label=\"{} (AUC = {:.3f})\".format(name, auc))\n\nplt.plot([0, 1], [0, 1], linestyle=\"--\", color=\"gray\")\nplt.title(\"ROC Curve Comparison\")\nplt.xlabel(\"False Positive Rate\")\nplt.ylabel(\"True Positive Rate\")\nplt.legend()\nplt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-08-26T11:29:30.181924Z","iopub.execute_input":"2026-08-26T11:29:30.182242Z","iopub.status.idle":"2026-08-26T11:29:30.362395Z","shell.execute_reply.started":"2026-08-26T11:29:30.182200Z","shell.execute_reply":"2026-08-26T11:29:30.361556Z"}},"outputs":[],"execution_count":null},{"id":"5636110a-7873-4ecc-89e2-091a81622580","cell_type":"markdown","source":"whichever model has the highest ROC-AUC and the curve closest to the top-left corner is the strongest baseline to carry forward into tuning","metadata":{}},{"id":"3f4f1f4d-86d9-48be-9d58-b964aac6088a","cell_type":"markdown","source":"---\n## 8. Hyperparameter Tuning with Optuna\n\nOptuna uses Bayesian optimization (TPE) instead of grid search: it builds a probabilistic model of which parameter regions score well and samples more from those regions on later trials, so it tends to find strong parameters in far fewer trials than a grid search would need. The best baseline model (XGBoost) is tuned here, optimizing ROC-AUC on a held-out validation split.","metadata":{}},{"id":"349ae437-c9de-400a-a310-c81f123c6973","cell_type":"code","source":"X_tr, X_val, y_tr, y_val = train_test_split(\n    X_train, y_train, test_size=0.2, random_state=SEED, stratify=y_train\n)\n\ndef objective(trial):\n    params = {\n        \"n_estimators\": trial.suggest_int(\"n_estimators\", 100, 600),\n        \"max_depth\": trial.suggest_int(\"max_depth\", 3, 9),\n        \"learning_rate\": trial.suggest_float(\"learning_rate\", 0.01, 0.3, log=True),\n        \"subsample\": trial.suggest_float(\"subsample\", 0.6, 1.0),\n        \"colsample_bytree\": trial.suggest_float(\"colsample_bytree\", 0.6, 1.0),\n        \"min_child_weight\": trial.suggest_int(\"min_child_weight\", 1, 10),\n        \"scale_pos_weight\": scale_pos_weight,\n        \"eval_metric\": \"logloss\",\n        \"random_state\": SEED\n    }\n    model = xgb.XGBClassifier(**params)\n    model.fit(X_tr, y_tr)\n    proba = model.predict_proba(X_val)[:, 1]\n    return roc_auc_score(y_val, proba)\n\nN_TRIALS = 40   # increase if you have time to spare\nstudy = optuna.create_study(direction=\"maximize\")\nstudy.optimize(objective, n_trials=N_TRIALS, show_progress_bar=False)\n\nprint(\"Best ROC-AUC on validation:\", round(study.best_value, 4))\nprint(\"Best params:\", study.best_params)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-08-26T11:29:30.363282Z","iopub.execute_input":"2026-08-26T11:29:30.363504Z","iopub.status.idle":"2026-08-26T11:29:37.503721Z","shell.execute_reply.started":"2026-08-26T11:29:30.363483Z","shell.execute_reply":"2026-08-26T11:29:37.502975Z"}},"outputs":[],"execution_count":null},{"id":"859a676c-4c87-48eb-8371-3ad3b967ddc6","cell_type":"code","source":"fig, ax = plt.subplots(figsize=(10, 5))\nvals = [t.value for t in study.trials]\nbest_so_far = [max(vals[:i + 1]) for i in range(len(vals))]\nax.plot(vals, color=\"steelblue\", alpha=0.5, lw=1, marker=\"o\", ms=3, label=\"Trial ROC-AUC\")\nax.plot(best_so_far, color=\"crimson\", lw=2, label=\"Best so far\")\nax.set_title(\"Optuna Search Progress\")\nax.set_xlabel(\"Trial\")\nax.set_ylabel(\"ROC-AUC\")\nax.legend()\nplt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-08-26T11:29:37.504885Z","iopub.execute_input":"2026-08-26T11:29:37.505614Z","iopub.status.idle":"2026-08-26T11:29:37.685450Z","shell.execute_reply.started":"2026-08-26T11:29:37.505553Z","shell.execute_reply":"2026-08-26T11:29:37.684380Z"}},"outputs":[],"execution_count":null},{"id":"26053ea1-fa3b-422f-8ed7-dbf581f706fb","cell_type":"code","source":"best_params = study.best_params\nbest_params.update({\"scale_pos_weight\": scale_pos_weight, \"eval_metric\": \"logloss\", \"random_state\": SEED})\n\ntuned_model = xgb.XGBClassifier(**best_params)\ntuned_model.fit(X_train, y_train)\n\nm, p = evaluate_model(\"XGBoost Tuned\", tuned_model, X_test, y_test)\nresults.append(m); proba_by_model[\"XGBoost Tuned\"] = p\n\nresults_df = pd.DataFrame(results).set_index(\"model\").round(4)\nresults_df","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-08-26T11:29:37.686695Z","iopub.execute_input":"2026-08-26T11:29:37.687291Z","iopub.status.idle":"2026-08-26T11:29:37.903549Z","shell.execute_reply.started":"2026-08-26T11:29:37.687253Z","shell.execute_reply":"2026-08-26T11:29:37.902826Z"}},"outputs":[],"execution_count":null},{"id":"8ea5de29-c737-4c5d-b4b2-6557e30797dc","cell_type":"markdown","source":"---\n## 9. Final Evaluation","metadata":{}},{"id":"35cad2c7-3c52-4c69-80a0-50b651491dec","cell_type":"code","source":"preds_tuned = tuned_model.predict(X_test)\ncm = confusion_matrix(y_test, preds_tuned)\n\nplt.figure(figsize=(6, 5))\nsns.heatmap(cm, annot=True, fmt=\"d\", cmap=\"Blues\",\n            xticklabels=[\"No Tsunami\", \"Tsunami\"], yticklabels=[\"No Tsunami\", \"Tsunami\"])\nplt.title(\"Confusion Matrix, Tuned XGBoost\")\nplt.xlabel(\"Predicted\")\nplt.ylabel(\"Actual\")\nplt.show()\n\nprint(classification_report(y_test, preds_tuned, target_names=[\"No Tsunami\", \"Tsunami\"]))","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-08-26T11:29:37.904487Z","iopub.execute_input":"2026-08-26T11:29:37.904797Z","iopub.status.idle":"2026-08-26T11:29:38.078390Z","shell.execute_reply.started":"2026-08-26T11:29:37.904773Z","shell.execute_reply":"2026-08-26T11:29:38.077629Z"}},"outputs":[],"execution_count":null},{"id":"5eca5d88-5ee7-4009-a691-fc288baf502e","cell_type":"code","source":"importances = pd.Series(tuned_model.feature_importances_, index=X.columns).sort_values(ascending=False)\n\nplt.figure(figsize=(9, 8))\nsns.barplot(x=importances.values[:15], y=importances.index[:15], color=\"steelblue\")\nplt.title(\"Top 15 Feature Importances, Tuned XGBoost\")\nplt.xlabel(\"Importance\")\nplt.ylabel(\"Feature\")\nplt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-08-26T11:29:38.079443Z","iopub.execute_input":"2026-08-26T11:29:38.079793Z","iopub.status.idle":"2026-08-26T11:29:38.293849Z","shell.execute_reply.started":"2026-08-26T11:29:38.079768Z","shell.execute_reply":"2026-08-26T11:29:38.292941Z"}},"outputs":[],"execution_count":null},{"id":"528bd5b5-fd4a-4e15-a856-6504edb377d8","cell_type":"markdown","source":"### Summary\n\nMagnitude and sig turned out to be the strongest signals for tsunami risk, both in the correlation analysis and in what the tuned model actually leaned on. PCA and K-Means show that this risk pattern shows up even without using the label at all, which is a good sign the signal is real and not something the model invented.\n\nThe tuned XGBoost model catches most tsunami events (recall 0.64, precision 0.75) but still misses a meaningful chunk, worth flagging since only 14 tsunami cases were in the test set, so this number isn't fully stable yet. Given the cost of missing a real tsunami, recall is the main thing to push higher going forward, whether through threshold tuning, more data, or added features like distance to known subduction zones.","metadata":{}}]}