{"cells":[{"cell_type":"code","execution_count":1,"id":"727c813e","metadata":{},"outputs":[],"source":"# Install required packages (safe to re-run; skips if already installed)\nimport subprocess, sys\npkgs = [\"xgboost\", \"lightgbm\", \"optuna\", \"shap\", \"scikit-learn\", \"matplotlib\", \"seaborn\"]\nfor pkg in pkgs:\n    try:\n        __import__(pkg.replace(\"-\", \"_\"))\n    except ImportError:\n        subprocess.check_call([sys.executable, \"-m\", \"pip\", \"install\", pkg, \"-q\"])\nprint(\"All packages ready.\")"},{"cell_type":"markdown","id":"7ebbbe0b","metadata":{},"source":"# Insurance Premium Prediction — Kaggle Playground S4E12\n\n**DSTA Final Project — Spring 2026**\n**Team:** Rafi · Maaz · Quaid\n\nA regression pipeline predicting `Premium Amount` from 19 customer/policy features.\n\n**Pipeline:**\n1. Environment setup & data loading\n2. Comprehensive EDA (12 charts) — **Rafi**\n3. Light feature engineering — **Maaz**\n4. Model-zoo — three self-contained individual solutions:\n    - Ridge regression baseline — **Rafi**\n    - LightGBM single — **Maaz**\n    - XGBoost single — **Quaid**\n5. Heavy feature engineering (target + count encoding, 15 COMBOs) — **Maaz**\n6. Optuna hyperparameter tuning — **Maaz**\n7. Final 5-fold XGB+LGB ensemble with blend search — **Quaid**\n8. Stability & diagnostics — **Quaid**\n9. SHAP explainability — **Rafi**\n10. Submission\n\n**Metric:** RMSE on `log1p(Premium Amount)`."},{"cell_type":"markdown","id":"0448c30f","metadata":{},"source":"## 1. Setup, Environment Detection & Data Loading"},{"cell_type":"code","execution_count":2,"id":"5b25206e","metadata":{},"outputs":[],"source":"import os\nimport gc\nimport time\nimport warnings\nimport subprocess\nfrom pathlib import Path\n\nimport numpy as np\nimport pandas as pd\nimport matplotlib.pyplot as plt\nimport seaborn as sns\n\nwarnings.filterwarnings(\"ignore\")\npd.set_option(\"display.max_columns\", 200)\npd.set_option(\"display.width\", 180)\nsns.set_style(\"whitegrid\")\n\nfrom sklearn.model_selection import KFold, train_test_split\nfrom sklearn.linear_model import Ridge\nfrom sklearn.preprocessing import StandardScaler\nfrom sklearn.metrics import mean_squared_error, mean_absolute_error, r2_score\nfrom sklearn.impute import SimpleImputer\n\nfrom xgboost import XGBRegressor\nimport lightgbm as lgb\n\ntry:\n    import optuna\n    optuna.logging.set_verbosity(optuna.logging.WARNING)\n    HAS_OPTUNA = True\nexcept Exception:\n    HAS_OPTUNA = False\n\ntry:\n    import shap\n    HAS_SHAP = True\nexcept Exception:\n    HAS_SHAP = False\n\ntry:\n    import cudf  # noqa\n    HAS_CUDF = True\nexcept Exception:\n    HAS_CUDF = False\n\ndef has_nvidia_gpu():\n    try:\n        subprocess.run([\"nvidia-smi\"], stdout=subprocess.PIPE, stderr=subprocess.PIPE, check=True)\n        return True\n    except Exception:\n        return False\n\nHAS_GPU = has_nvidia_gpu()\nRANDOM_STATE = 42\nnp.random.seed(RANDOM_STATE)\n\nCHARTS_DIR = Path(\"notebook_charts_insurance\")\nCHARTS_DIR.mkdir(exist_ok=True)\n\ndef save_fig(name, dpi=120):\n    plt.tight_layout()\n    plt.savefig(CHARTS_DIR / f\"{name}.png\", dpi=dpi, bbox_inches=\"tight\")\n    plt.show()\n\nprint(f\"HAS_GPU    = {HAS_GPU}\")\nprint(f\"HAS_CUDF   = {HAS_CUDF}\")\nprint(f\"HAS_OPTUNA = {HAS_OPTUNA}\")\nprint(f\"HAS_SHAP   = {HAS_SHAP}\")"},{"cell_type":"code","execution_count":3,"metadata":{},"outputs":[],"source":"# Data loading — Kaggle path first, else playground-series-s4e12/, else project root.\nCANDIDATE_DIRS = [\n    Path(\"/kaggle/input/competitions/playground-series-s4e12\"),\n    Path(\"playground-series-s4e12\"),\n    Path(\".\"),\n]\nTRAIN_PATH, TEST_PATH = None, None\nfor d in CANDIDATE_DIRS:\n    if (d / \"train.csv\").exists() and (d / \"test.csv\").exists():\n        TRAIN_PATH = d / \"train.csv\"\n        TEST_PATH = d / \"test.csv\"\n        break\nassert TRAIN_PATH is not None, \"Could not find train.csv / test.csv\"\nprint(f\"Loading from: {TRAIN_PATH.parent}\")\n\ntrain_raw = pd.read_csv(TRAIN_PATH)\ntest_raw = pd.read_csv(TEST_PATH)\nprint(f\"train shape: {train_raw.shape}\")\nprint(f\"test  shape: {test_raw.shape}\")\ntrain_raw.head()"},{"cell_type":"markdown","metadata":{},"source":"## 2. Exploratory Data Analysis — *Rafi*\n\nTwelve charts covering the target distribution, missingness, numeric and\ncategorical structure, temporal behaviour, and interactions. The goal is to\nsurface every signal we will later exploit in feature engineering."},{"cell_type":"code","execution_count":4,"metadata":{},"outputs":[],"source":"TARGET = \"Premium Amount\"\nprint(f\"Target stats:\")\nprint(train_raw[TARGET].describe())\nprint(f\"\\nSkew (raw):   {train_raw[TARGET].skew():.3f}\")\nprint(f\"Skew (log1p): {np.log1p(train_raw[TARGET]).skew():.3f}\")"},{"cell_type":"code","execution_count":null,"metadata":{},"outputs":[],"source":"# 2.1 Target distribution — raw vs log1p\nfig, axes = plt.subplots(1, 2, figsize=(14, 4.5))\naxes[0].hist(train_raw[TARGET], bins=60, color=\"#1E2761\", alpha=0.85)\naxes[0].set_title(\"Raw Premium Amount\", fontweight=\"bold\")\naxes[0].set_xlabel(\"Premium Amount\"); axes[0].set_ylabel(\"Count\")\naxes[1].hist(np.log1p(train_raw[TARGET]), bins=60, color=\"#F96167\", alpha=0.85)\naxes[1].set_title(\"log1p(Premium Amount) — near-normal\", fontweight=\"bold\")\naxes[1].set_xlabel(\"log1p(Premium)\")\nplt.suptitle(\"Target Distribution — justifies log-transform\", fontsize=14, fontweight=\"bold\")\nsave_fig(\"01_target_distribution\")"},{"cell_type":"markdown","metadata":{},"source":"**Takeaway:** The raw target is heavily right-skewed (skew ≈ 3.4). Applying `log1p` pulls it close to normal — we train on the log scale and exponentiate back for submission."},{"cell_type":"code","execution_count":null,"metadata":{},"outputs":[],"source":"# 2.2 Missingness bar chart\nmiss = train_raw.isnull().sum().sort_values(ascending=False)\nmiss = miss[miss > 0]\nplt.figure(figsize=(10, 4.5))\nmiss.plot(kind=\"bar\", color=\"#F96167\")\nplt.title(\"Missing values by column (train)\", fontweight=\"bold\")\nplt.ylabel(\"# missing rows\")\nplt.xticks(rotation=45, ha=\"right\")\nsave_fig(\"02_missingness\")\nprint(miss)"},{"cell_type":"markdown","metadata":{},"source":"**Takeaway:** Occupation is missing ~30%. Credit Score, Vehicle Age, Age, Annual Income, Health Score, Previous Claims all have meaningful missingness — we add boolean missing-indicator flags for the top columns."},{"cell_type":"code","execution_count":null,"metadata":{},"outputs":[],"source":"# 2.3 Numeric feature distributions\nNUMERIC_COLS = [\"Age\", \"Annual Income\", \"Health Score\", \"Credit Score\",\n                \"Vehicle Age\", \"Insurance Duration\", \"Previous Claims\",\n                \"Number of Dependents\"]\nfig, axes = plt.subplots(2, 4, figsize=(16, 7))\nfor ax, col in zip(axes.ravel(), NUMERIC_COLS):\n    data = train_raw[col].dropna()\n    ax.hist(data, bins=40, color=\"#028090\", alpha=0.85)\n    ax.set_title(col, fontweight=\"bold\", fontsize=10)\nplt.suptitle(\"Numeric feature distributions\", fontsize=14, fontweight=\"bold\")\nsave_fig(\"03_numeric_distributions\")"},{"cell_type":"code","execution_count":null,"metadata":{},"outputs":[],"source":"# 2.4 Numeric correlation with log-target\ny_log = np.log1p(train_raw[TARGET])\ncorrs = {}\nfor col in NUMERIC_COLS:\n    corrs[col] = train_raw[col].fillna(train_raw[col].median()).corr(y_log)\ncorr_s = pd.Series(corrs).sort_values()\nplt.figure(figsize=(9, 4.5))\ncolors = [\"#F96167\" if v < 0 else \"#028090\" for v in corr_s.values]\nplt.barh(corr_s.index, corr_s.values, color=colors)\nplt.title(\"Pearson correlation of numeric features with log1p(Premium)\", fontweight=\"bold\")\nplt.axvline(0, color=\"k\", linewidth=0.8)\nsave_fig(\"04_numeric_correlation\")\nprint(corr_s)"},{"cell_type":"markdown","metadata":{},"source":"**Takeaway:** Numeric correlations are weak individually (|r| < 0.05). The signal lives in interactions and categorical structure — we will engineer ratios and lean on gradient-boosted trees."},{"cell_type":"code","execution_count":null,"metadata":{},"outputs":[],"source":"# 2.5 Categorical cardinality + mean target by category\nCAT_COLS = [\"Gender\", \"Marital Status\", \"Education Level\", \"Occupation\",\n            \"Location\", \"Policy Type\", \"Customer Feedback\", \"Smoking Status\",\n            \"Exercise Frequency\", \"Property Type\"]\nfig, axes = plt.subplots(3, 4, figsize=(16, 11))\naxes = axes.ravel()\nfor i, col in enumerate(CAT_COLS):\n    ax = axes[i]\n    grp = train_raw.groupby(col, dropna=False)[TARGET].mean().sort_values()\n    grp.plot(kind=\"bar\", ax=ax, color=\"#1E2761\")\n    ax.set_title(f\"{col}  (n={train_raw[col].nunique(dropna=False)})\", fontweight=\"bold\", fontsize=10)\n    ax.set_ylabel(\"mean Premium\")\n    ax.tick_params(axis=\"x\", rotation=30, labelsize=8)\nfor j in range(len(CAT_COLS), len(axes)):\n    axes[j].axis(\"off\")\nplt.suptitle(\"Mean Premium by categorical value\", fontsize=14, fontweight=\"bold\")\nsave_fig(\"05_categorical_mean_target\")"},{"cell_type":"markdown","metadata":{},"source":"**Takeaway:** Customer Feedback, Policy Type and Smoking Status show the clearest mean-premium spreads — these will drive early splits in the gradient-boosted models."},{"cell_type":"code","execution_count":null,"metadata":{},"outputs":[],"source":"# 2.6 Policy Start Date — monthly volume + mean premium\npdate = pd.to_datetime(train_raw[\"Policy Start Date\"])\ntmp = pd.DataFrame({\"ym\": pdate.dt.to_period(\"M\").astype(str), \"premium\": train_raw[TARGET]})\nmonthly = tmp.groupby(\"ym\")[\"premium\"].agg([\"count\", \"mean\"]).reset_index()\nfig, ax1 = plt.subplots(figsize=(12, 4.5))\nax1.bar(monthly[\"ym\"], monthly[\"count\"], color=\"#CADCFC\", label=\"volume\")\nax1.set_ylabel(\"# policies\", color=\"#1E2761\")\nax1.tick_params(axis=\"x\", rotation=90, labelsize=7)\nax2 = ax1.twinx()\nax2.plot(monthly[\"ym\"], monthly[\"mean\"], color=\"#F96167\", linewidth=2, label=\"mean premium\")\nax2.set_ylabel(\"mean Premium\", color=\"#F96167\")\nplt.title(\"Policy Start Date — volume vs mean premium by month\", fontweight=\"bold\")\nsave_fig(\"06_policy_date_trend\")"},{"cell_type":"code","execution_count":null,"metadata":{},"outputs":[],"source":"# 2.7 Age-bucket mean premium\nage_bins = [0, 25, 35, 45, 55, 65, 100]\nage_lbl = [\"<25\", \"25-34\", \"35-44\", \"45-54\", \"55-64\", \"65+\"]\nab = pd.cut(train_raw[\"Age\"], bins=age_bins, labels=age_lbl, right=False)\nag = train_raw.groupby(ab, observed=True)[TARGET].mean()\nplt.figure(figsize=(9, 4))\nag.plot(kind=\"bar\", color=\"#1E2761\")\nplt.title(\"Mean Premium by Age bucket\", fontweight=\"bold\")\nplt.ylabel(\"mean Premium\")\nplt.xticks(rotation=0)\nsave_fig(\"07_age_bucket\")"},{"cell_type":"code","execution_count":null,"metadata":{},"outputs":[],"source":"# 2.8 Income-decile mean premium\ninc = train_raw[\"Annual Income\"].fillna(train_raw[\"Annual Income\"].median())\ndec = pd.qcut(inc, 10, labels=False, duplicates=\"drop\")\nidf = pd.DataFrame({\"decile\": dec, \"p\": train_raw[TARGET]}).groupby(\"decile\")[\"p\"].mean()\nplt.figure(figsize=(9, 4))\nidf.plot(marker=\"o\", color=\"#028090\", linewidth=2)\nplt.title(\"Mean Premium by Annual-Income decile\", fontweight=\"bold\")\nplt.ylabel(\"mean Premium\"); plt.xlabel(\"income decile\")\nsave_fig(\"08_income_decile\")"},{"cell_type":"code","execution_count":null,"metadata":{},"outputs":[],"source":"# 2.9 Health Score vs mean premium\nhs = train_raw[\"Health Score\"].fillna(train_raw[\"Health Score\"].median())\nhb = pd.cut(hs, bins=20)\nhdf = train_raw.groupby(hb, observed=True)[TARGET].mean()\nplt.figure(figsize=(10, 4))\nhdf.plot(color=\"#F96167\", linewidth=2, marker=\"o\", markersize=4)\nplt.title(\"Mean Premium by Health Score bucket\", fontweight=\"bold\")\nplt.ylabel(\"mean Premium\"); plt.xlabel(\"Health Score bucket\")\nplt.xticks(rotation=45, ha=\"right\", fontsize=7)\nsave_fig(\"09_health_score\")"},{"cell_type":"code","execution_count":null,"metadata":{},"outputs":[],"source":"# 2.10 Previous Claims (capped 3+) vs mean premium\npc = train_raw[\"Previous Claims\"].fillna(0).clip(upper=3).astype(int)\npdf = pd.DataFrame({\"pc\": pc, \"p\": train_raw[TARGET]}).groupby(\"pc\")[\"p\"].agg([\"mean\", \"count\"])\nfig, ax1 = plt.subplots(figsize=(8, 4))\nax1.bar(pdf.index, pdf[\"mean\"], color=\"#1E2761\")\nax1.set_ylabel(\"mean Premium\", color=\"#1E2761\")\nax1.set_xlabel(\"Previous Claims (capped at 3+)\")\nax2 = ax1.twinx()\nax2.plot(pdf.index, pdf[\"count\"], color=\"#F96167\", linewidth=2, marker=\"o\")\nax2.set_ylabel(\"count\", color=\"#F96167\")\nplt.title(\"Previous Claims vs mean Premium\", fontweight=\"bold\")\nsave_fig(\"10_prev_claims\")"},{"cell_type":"code","execution_count":null,"metadata":{},"outputs":[],"source":"# 2.11 Smoking × Exercise heatmap of mean premium\npiv = train_raw.pivot_table(index=\"Smoking Status\", columns=\"Exercise Frequency\",\n                            values=TARGET, aggfunc=\"mean\")\nplt.figure(figsize=(8, 4))\nsns.heatmap(piv, annot=True, fmt=\".0f\", cmap=\"RdBu_r\", center=train_raw[TARGET].mean())\nplt.title(\"Mean Premium: Smoking × Exercise Frequency\", fontweight=\"bold\")\nsave_fig(\"11_smoking_exercise_heatmap\")"},{"cell_type":"code","execution_count":null,"metadata":{},"outputs":[],"source":"# 2.12 Correlation heatmap of numerics\nnum_df = train_raw[NUMERIC_COLS + [TARGET]].copy()\nnum_df[TARGET] = np.log1p(num_df[TARGET])\nplt.figure(figsize=(9, 7))\nsns.heatmap(num_df.corr(), annot=True, fmt=\".2f\", cmap=\"RdBu_r\", center=0,\n            vmin=-0.1, vmax=0.1)\nplt.title(\"Numeric correlation heatmap (target = log1p)\", fontweight=\"bold\")\nsave_fig(\"12_correlation_heatmap\")"},{"cell_type":"markdown","metadata":{},"source":"### EDA summary (Rafi)\n\n- Target is heavily right-skewed → train on `log1p`.\n- Missingness is concentrated in `Occupation` (30%), with moderate gaps in Credit Score, Vehicle Age, Age, Income, Health Score, Previous Claims → add missing-flag features.\n- Numeric Pearson correlation with the log target is tiny (|r| < 0.05) individually. No single dominant driver — this is why aggressive feature engineering (interactions, target encoding) pays off.\n- The biggest categorical signals are Customer Feedback, Policy Type, Smoking Status. Smoker × Exercise-Frequency interaction is clearly non-additive.\n- Policy Start Date shows mild seasonal variation. Age and Income relationships with premium are non-linear — gradient-boosted trees will capture this without manual splines."},{"cell_type":"markdown","metadata":{},"source":"## 3. Light Feature Engineering — *Maaz*\n\nDate decomposition, cyclic encoding, ratio features, and missing-indicator\nflags. This \"light\" feature set feeds the three member baselines in\nsection 4; the heavy target-encoding pipeline is built later in section 5."},{"cell_type":"code","execution_count":null,"metadata":{},"outputs":[],"source":"def engineer_light(df):\n    df = df.copy()\n    pd_dt = pd.to_datetime(df[\"Policy Start Date\"])\n    df[\"year\"] = pd_dt.dt.year.astype(\"float32\")\n    df[\"month\"] = pd_dt.dt.month.astype(\"float32\")\n    df[\"day\"] = pd_dt.dt.day.astype(\"float32\")\n    df[\"dow\"] = pd_dt.dt.dayofweek.astype(\"float32\")\n    df[\"woy\"] = pd_dt.dt.isocalendar().week.astype(\"int32\").astype(\"float32\")\n    df[\"seconds\"] = (pd_dt.astype(\"int64\") // 10**9).astype(\"float32\")\n    df[\"month_sin\"] = np.sin(2 * np.pi * df[\"month\"] / 12).astype(\"float32\")\n    df[\"month_cos\"] = np.cos(2 * np.pi * df[\"month\"] / 12).astype(\"float32\")\n    df[\"quarter\"] = ((df[\"month\"] - 1) // 3 + 1).astype(\"float32\")\n    df[\"is_weekend_start\"] = (df[\"dow\"] >= 5).astype(\"float32\")\n    # Missing-indicator flags\n    for c in [\"Occupation\", \"Credit Score\", \"Vehicle Age\", \"Age\",\n              \"Annual Income\", \"Health Score\", \"Previous Claims\"]:\n        df[f\"missing_{c.replace(' ', '_')}\"] = df[c].isna().astype(\"float32\")\n    # Ratios (Laplace-smoothed where needed)\n    df[\"income_per_age\"] = (df[\"Annual Income\"].fillna(df[\"Annual Income\"].median())\n                            / (df[\"Age\"].fillna(df[\"Age\"].median()) + 1)).astype(\"float32\")\n    df[\"health_income_ratio\"] = (df[\"Health Score\"].fillna(50)\n                                  / np.log1p(df[\"Annual Income\"].fillna(df[\"Annual Income\"].median()))\n                                 ).astype(\"float32\")\n    df[\"claims_rate\"] = (df[\"Previous Claims\"].fillna(0)\n                          / (df[\"Insurance Duration\"].fillna(1) + 1)).astype(\"float32\")\n    df[\"has_prev_claim\"] = (df[\"Previous Claims\"].fillna(0) > 0).astype(\"float32\")\n    df[\"credit_income_ratio\"] = (df[\"Credit Score\"].fillna(train_raw[\"Credit Score\"].median())\n                                  / np.log1p(df[\"Annual Income\"].fillna(df[\"Annual Income\"].median()))\n                                 ).astype(\"float32\")\n    df[\"vehicle_age_ratio\"] = (df[\"Vehicle Age\"].fillna(0) / (df[\"Age\"].fillna(df[\"Age\"].median()) + 1)\n                               ).astype(\"float32\")\n    df[\"age_bucket\"] = pd.cut(df[\"Age\"].fillna(df[\"Age\"].median()),\n                              bins=[0, 25, 35, 45, 55, 65, 100], labels=False).astype(\"float32\")\n    return df\n\ntrain_fe = engineer_light(train_raw)\ntest_fe = engineer_light(test_raw)\nprint(\"After light FE — train:\", train_fe.shape, \"test:\", test_fe.shape)"},{"cell_type":"code","execution_count":null,"metadata":{},"outputs":[],"source":"# Encode categoricals with factorize, impute numerics with median.\ndef prep_light(tr, te):\n    tr = tr.copy(); te = te.copy()\n    tr[\"y\"] = np.log1p(tr[TARGET]).astype(\"float32\")\n    combined = pd.concat([tr.drop(columns=[\"y\", TARGET]), te], axis=0, ignore_index=True)\n    # Drop non-feature cols\n    DROP = [\"id\", \"Policy Start Date\", TARGET]\n    feats = [c for c in combined.columns if c not in DROP]\n    for c in feats:\n        if combined[c].dtype == \"object\":\n            combined[c] = combined[c].fillna(\"NAN\").astype(str)\n            codes, _ = pd.factorize(combined[c], sort=True)\n            combined[c] = codes.astype(\"int32\")\n        elif combined[c].dtype == \"float64\":\n            combined[c] = combined[c].astype(\"float32\")\n    # Impute remaining NaNs with column median\n    for c in feats:\n        if combined[c].isna().any():\n            combined[c] = combined[c].fillna(combined[c].median())\n    n_tr = len(tr)\n    tr_out = pd.concat([tr[[\"id\", \"y\"]].reset_index(drop=True),\n                        combined.iloc[:n_tr][feats].reset_index(drop=True)], axis=1)\n    te_out = pd.concat([te[[\"id\"]].reset_index(drop=True),\n                        combined.iloc[n_tr:][feats].reset_index(drop=True)], axis=1)\n    return tr_out, te_out, feats\n\ntrain_L, test_L, LIGHT_FEATURES = prep_light(train_fe, test_fe)\nprint(f\"Light feature count: {len(LIGHT_FEATURES)}\")\nprint(LIGHT_FEATURES)"},{"cell_type":"markdown","metadata":{},"source":"## 4. Model Zoo — Three Self-Contained Individual Solutions\n\nEach team member developed a standalone baseline. These prove every model in\nthe final ensemble is our own work (self-contained), and they quantify how\nmuch the ensemble actually gains over a single model.\n\n- **Ridge regression** — *Rafi* (interpretable linear baseline).\n- **LightGBM single** — *Maaz*.\n- **XGBoost single** — *Quaid*.\n\nCV: common 5-fold KFold split on the full 1.2M rows, RMSE on the log target."},{"cell_type":"code","execution_count":null,"metadata":{},"outputs":[],"source":"FOLDS = 5\nkf = KFold(n_splits=FOLDS, shuffle=True, random_state=RANDOM_STATE)\nX_light = train_L[LIGHT_FEATURES].values.astype(np.float32)\ny_light = train_L[\"y\"].values.astype(np.float32)\nXt_light = test_L[LIGHT_FEATURES].values.astype(np.float32)\n\ndef rmse(a, b):\n    return float(np.sqrt(mean_squared_error(a, b)))\n\nzoo_results = {}"},{"cell_type":"code","execution_count":null,"metadata":{},"outputs":[],"source":"# 4.1 Ridge — Rafi\nprint(\"=== Ridge (Rafi) ===\")\nt0 = time.time()\noof_ridge = np.zeros(len(train_L), dtype=np.float32)\npred_ridge = np.zeros(len(test_L), dtype=np.float32)\nscaler = StandardScaler()\nX_light_s = scaler.fit_transform(X_light)\nXt_light_s = scaler.transform(Xt_light)\nfor fi, (tr, va) in enumerate(kf.split(X_light_s)):\n    m = Ridge(alpha=1.0, random_state=RANDOM_STATE)\n    m.fit(X_light_s[tr], y_light[tr])\n    oof_ridge[va] = m.predict(X_light_s[va])\n    pred_ridge += m.predict(Xt_light_s) / FOLDS\nridge_rmse = rmse(y_light, oof_ridge)\nridge_time = time.time() - t0\nzoo_results[\"Ridge (Rafi)\"] = dict(rmse=ridge_rmse, mae=mean_absolute_error(y_light, oof_ridge),\n                                    r2=r2_score(y_light, oof_ridge), time=ridge_time)\nprint(f\"Ridge  OOF RMSE = {ridge_rmse:.5f}   time = {ridge_time:.1f}s\")"},{"cell_type":"code","execution_count":null,"metadata":{},"outputs":[],"source":"# 4.2 LightGBM — Quaid\nprint(\"=== LightGBM single (Quaid) ===\")\nt0 = time.time()\noof_lgb_single = np.zeros(len(train_L), dtype=np.float32)\npred_lgb_single = np.zeros(len(test_L), dtype=np.float32)\nlgb_base_params = dict(objective=\"regression\", metric=\"rmse\", verbosity=-1,\n                       n_estimators=800, learning_rate=0.05, num_leaves=128,\n                       min_child_samples=25, subsample=0.85, subsample_freq=1,\n                       colsample_bytree=0.85, reg_alpha=0.05, reg_lambda=0.5,\n                       random_state=RANDOM_STATE, n_jobs=-1)\nfor fi, (tr, va) in enumerate(kf.split(X_light)):\n    m = lgb.LGBMRegressor(**lgb_base_params)\n    m.fit(X_light[tr], y_light[tr], eval_set=[(X_light[va], y_light[va])],\n          callbacks=[lgb.early_stopping(30, verbose=False)])\n    oof_lgb_single[va] = m.predict(X_light[va])\n    pred_lgb_single += m.predict(Xt_light) / FOLDS\nlgb_single_rmse = rmse(y_light, oof_lgb_single)\nlgb_single_time = time.time() - t0\nzoo_results[\"LightGBM (Maaz)\"] = dict(rmse=lgb_single_rmse, mae=mean_absolute_error(y_light, oof_lgb_single),\n                                       r2=r2_score(y_light, oof_lgb_single), time=lgb_single_time)\nprint(f\"LGB    OOF RMSE = {lgb_single_rmse:.5f}   time = {lgb_single_time:.1f}s\")"},{"cell_type":"code","execution_count":null,"metadata":{},"outputs":[],"source":"# 4.3 XGBoost — Rafi\nprint(\"=== XGBoost single (Rafi) ===\")\nt0 = time.time()\noof_xgb_single = np.zeros(len(train_L), dtype=np.float32)\npred_xgb_single = np.zeros(len(test_L), dtype=np.float32)\nxgb_base_params = dict(max_depth=7, colsample_bytree=0.85, subsample=0.85,\n                        n_estimators=800, learning_rate=0.05,\n                        early_stopping_rounds=30, eval_metric=\"rmse\",\n                        random_state=RANDOM_STATE, tree_method=\"hist\",\n                        device=\"cuda\" if HAS_GPU else \"cpu\")\nfor fi, (tr, va) in enumerate(kf.split(X_light)):\n    m = XGBRegressor(**xgb_base_params)\n    m.fit(X_light[tr], y_light[tr], eval_set=[(X_light[va], y_light[va])], verbose=False)\n    oof_xgb_single[va] = m.predict(X_light[va])\n    pred_xgb_single += m.predict(Xt_light) / FOLDS\nxgb_single_rmse = rmse(y_light, oof_xgb_single)\nxgb_single_time = time.time() - t0\nzoo_results[\"XGBoost (Quaid)\"] = dict(rmse=xgb_single_rmse, mae=mean_absolute_error(y_light, oof_xgb_single),\n                                       r2=r2_score(y_light, oof_xgb_single), time=xgb_single_time)\nprint(f\"XGB    OOF RMSE = {xgb_single_rmse:.5f}   time = {xgb_single_time:.1f}s\")"},{"cell_type":"code","execution_count":null,"metadata":{},"outputs":[],"source":"# 4.4 Comparison table + chart\nzoo_df = pd.DataFrame(zoo_results).T\nzoo_df[\"rmse\"] = zoo_df[\"rmse\"].astype(float)\nzoo_df[\"mae\"] = zoo_df[\"mae\"].astype(float)\nzoo_df[\"r2\"] = zoo_df[\"r2\"].astype(float)\nzoo_df[\"time\"] = zoo_df[\"time\"].astype(float)\ndisplay(zoo_df.style.format({\"rmse\": \"{:.5f}\", \"mae\": \"{:.5f}\", \"r2\": \"{:.4f}\", \"time\": \"{:.1f}s\"}))\n\nfig, axes = plt.subplots(1, 2, figsize=(13, 4.2))\nzoo_df[\"rmse\"].plot(kind=\"bar\", ax=axes[0], color=[\"#1E2761\", \"#F96167\", \"#028090\"])\naxes[0].set_title(\"OOF RMSE by model\", fontweight=\"bold\"); axes[0].set_ylabel(\"RMSE\")\naxes[0].tick_params(axis=\"x\", rotation=15)\nzoo_df[\"time\"].plot(kind=\"bar\", ax=axes[1], color=[\"#1E2761\", \"#F96167\", \"#028090\"])\naxes[1].set_title(\"Training time (s)\", fontweight=\"bold\")\naxes[1].tick_params(axis=\"x\", rotation=15)\nsave_fig(\"13_model_zoo_comparison\")"},{"cell_type":"markdown","metadata":{},"source":"**Takeaway:** On light features the two tree-based models (LGB / XGB) clearly beat the Ridge baseline. XGB and LGB are near-tied — a textbook case for a blended ensemble. The heavy FE + Optuna-tuned ensemble in section 7 will push the RMSE further."},{"cell_type":"markdown","metadata":{},"source":"## 5. Heavy Feature Engineering — *Maaz*\n\nPer-fold target encoding (mean, median, min, max, std, nunique) and count\nencoding over **15 domain-motivated COMBOs** plus every single raw feature.\nEncoding is computed **inside each CV fold** with a 5-way inner KFold to\nprevent target leakage. Total engineered features: ~**150**."},{"cell_type":"code","execution_count":null,"metadata":{},"outputs":[],"source":"def gpu_groupby_agg(pdf, cols, target, agg):\n    use = cols + [target]\n    if agg == \"nunique\":\n        s = pdf.groupby(cols, dropna=False)[target].agg([\"nunique\", \"count\"]).reset_index()\n    elif agg == \"std\":\n        s = pdf.groupby(cols, dropna=False)[target].agg([\"std\", \"count\"]).reset_index()\n    else:\n        s = pdf.groupby(cols, dropna=False)[target].agg([agg, \"count\"]).reset_index()\n    s.columns = cols + [agg, \"count\"]\n    return s\n\ndef target_encode(train_df, valid_df, test_df, cols, target=\"y\", kfold=5, smooth=30, agg=\"mean\"):\n    feat = \"TE_\" + agg.upper() + \"_\" + \"_\".join(c.replace(\" \", \"_\") for c in cols)\n    train_df = train_df.copy(); valid_df = valid_df.copy(); test_df = test_df.copy()\n    train_df[\"_k\"] = np.arange(len(train_df)) % kfold\n    train_df[feat] = 0.0\n    if agg == \"mean\":  gs = train_df[target].mean()\n    elif agg == \"median\": gs = train_df[target].median()\n    elif agg == \"min\": gs = train_df[target].min()\n    elif agg == \"max\": gs = train_df[target].max()\n    elif agg == \"std\": gs = float(train_df[target].std())\n    elif agg == \"nunique\": gs = 0.0\n    for f in range(kfold):\n        fit = train_df[train_df[\"_k\"] != f]\n        mask = train_df[\"_k\"] == f\n        stats = gpu_groupby_agg(fit, cols, target, agg)\n        if agg in (\"nunique\", \"std\"):\n            stats[\"v\"] = stats[agg] / (stats[\"count\"] + 1.0)\n        else:\n            stats[\"v\"] = (stats[agg] * stats[\"count\"] + gs * smooth) / (stats[\"count\"] + smooth)\n        enc = train_df.loc[mask, cols].merge(stats[cols + [\"v\"]], on=cols, how=\"left\")[\"v\"].fillna(gs).values\n        train_df.loc[mask, feat] = enc\n    full = gpu_groupby_agg(train_df, cols, target, agg)\n    if agg in (\"nunique\", \"std\"):\n        full[\"v\"] = full[agg] / (full[\"count\"] + 1.0)\n    else:\n        full[\"v\"] = (full[agg] * full[\"count\"] + gs * smooth) / (full[\"count\"] + smooth)\n    for name, d in [(\"v\", valid_df), (\"t\", test_df)]:\n        e = d[cols].merge(full[cols + [\"v\"]], on=cols, how=\"left\")[\"v\"].fillna(gs).astype(\"float32\").values\n        (valid_df if name == \"v\" else test_df)[feat] = e\n    train_df[feat] = train_df[feat].astype(\"float32\")\n    train_df.drop(columns=[\"_k\"], inplace=True)\n    return train_df, valid_df, test_df\n\ndef count_encode(train_df, valid_df, test_df, cols):\n    feat = \"CE_\" + \"_\".join(c.replace(\" \", \"_\") for c in cols)\n    cnt = train_df.groupby(cols, dropna=False).size().reset_index(name=feat)\n    out = []\n    for d in [train_df, valid_df, test_df]:\n        m = d.merge(cnt, on=cols, how=\"left\")\n        m[feat] = m[feat].fillna(0).astype(\"int32\")\n        out.append(m)\n    return out[0], out[1], out[2]"},{"cell_type":"code","execution_count":null,"metadata":{},"outputs":[],"source":"# 15 domain-motivated COMBOS (trimmed from Original's 30 for speed)\nCOMBOS = [\n    [\"Annual Income\", \"Health Score\"],\n    [\"Credit Score\", \"Health Score\"],\n    [\"Customer Feedback\", \"Smoking Status\"],\n    [\"Health Score\", \"Smoking Status\"],\n    [\"Age\", \"Health Score\"],\n    [\"Exercise Frequency\", \"Smoking Status\"],\n    [\"Health Score\", \"Occupation\"],\n    [\"Policy Type\", \"Location\"],\n    [\"Credit Score\", \"Annual Income\"],\n    [\"Age\", \"Smoking Status\"],\n    [\"Health Score\", \"Vehicle Age\"],\n    [\"Age\", \"Credit Score\"],\n    [\"Education Level\", \"Marital Status\"],\n    [\"Previous Claims\", \"Vehicle Age\"],\n    [\"Health Score\", \"Insurance Duration\"],\n]\nprint(f\"COMBOS: {len(COMBOS)}\")"},{"cell_type":"markdown","metadata":{},"source":"## 6. Optuna Hyperparameter Tuning — *Maaz*\n\n20 trials each for LightGBM and XGBoost, TPE sampler, 3-fold CV on a 200k\nstratified sample to keep the search tractable. Objective: minimise CV RMSE\non the light-feature representation (heavy FE inherits these params)."},{"cell_type":"code","execution_count":null,"metadata":{},"outputs":[],"source":"# 200k stratified sample by log-target decile\nfrom sklearn.model_selection import StratifiedKFold\nsample_n = 200_000\ndec = pd.qcut(y_light, 10, labels=False, duplicates=\"drop\")\nidx = []\nfor d in np.unique(dec):\n    ii = np.where(dec == d)[0]\n    np.random.shuffle(ii)\n    idx.append(ii[:sample_n // 10])\nidx = np.concatenate(idx)\nX_opt = X_light[idx]; y_opt = y_light[idx]\nprint(f\"Optuna sample: {X_opt.shape}\")"},{"cell_type":"code","execution_count":null,"metadata":{},"outputs":[],"source":"def objective_lgb(trial):\n    params = dict(\n        objective=\"regression\", metric=\"rmse\", verbosity=-1,\n        num_leaves=trial.suggest_int(\"num_leaves\", 31, 256),\n        learning_rate=trial.suggest_float(\"learning_rate\", 5e-3, 1e-1, 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        reg_alpha=trial.suggest_float(\"reg_alpha\", 1e-4, 1.0, log=True),\n        reg_lambda=trial.suggest_float(\"reg_lambda\", 1e-4, 1.0, log=True),\n        min_child_samples=trial.suggest_int(\"min_child_samples\", 10, 100),\n        n_estimators=600, subsample_freq=1,\n        random_state=RANDOM_STATE, n_jobs=-1,\n    )\n    kfo = KFold(n_splits=3, shuffle=True, random_state=RANDOM_STATE)\n    scores = []\n    for tr, va in kfo.split(X_opt):\n        m = lgb.LGBMRegressor(**params)\n        m.fit(X_opt[tr], y_opt[tr], eval_set=[(X_opt[va], y_opt[va])],\n              callbacks=[lgb.early_stopping(25, verbose=False)])\n        scores.append(rmse(y_opt[va], m.predict(X_opt[va])))\n    return float(np.mean(scores))\n\nif HAS_OPTUNA:\n    t0 = time.time()\n    study_lgb = optuna.create_study(direction=\"minimize\", sampler=optuna.samplers.TPESampler(seed=RANDOM_STATE))\n    study_lgb.optimize(objective_lgb, n_trials=20, show_progress_bar=False)\n    print(f\"LGB Optuna: best RMSE = {study_lgb.best_value:.5f}   time = {time.time()-t0:.1f}s\")\n    print(study_lgb.best_params)\n    best_lgb_params = study_lgb.best_params\nelse:\n    best_lgb_params = dict(num_leaves=128, learning_rate=0.03, subsample=0.85,\n                           colsample_bytree=0.85, reg_alpha=0.05, reg_lambda=0.5,\n                           min_child_samples=25)"},{"cell_type":"code","execution_count":null,"metadata":{},"outputs":[],"source":"def objective_xgb(trial):\n    params = dict(\n        eval_metric=\"rmse\", tree_method=\"hist\",\n        device=\"cuda\" if HAS_GPU else \"cpu\",\n        max_depth=trial.suggest_int(\"max_depth\", 4, 12),\n        learning_rate=trial.suggest_float(\"learning_rate\", 5e-3, 1e-1, 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, 20),\n        reg_alpha=trial.suggest_float(\"reg_alpha\", 1e-4, 1.0, log=True),\n        reg_lambda=trial.suggest_float(\"reg_lambda\", 1e-4, 1.0, log=True),\n        n_estimators=600, early_stopping_rounds=25,\n        random_state=RANDOM_STATE,\n    )\n    kfo = KFold(n_splits=3, shuffle=True, random_state=RANDOM_STATE)\n    scores = []\n    for tr, va in kfo.split(X_opt):\n        m = XGBRegressor(**params)\n        m.fit(X_opt[tr], y_opt[tr], eval_set=[(X_opt[va], y_opt[va])], verbose=False)\n        scores.append(rmse(y_opt[va], m.predict(X_opt[va])))\n    return float(np.mean(scores))\n\nif HAS_OPTUNA:\n    t0 = time.time()\n    study_xgb = optuna.create_study(direction=\"minimize\", sampler=optuna.samplers.TPESampler(seed=RANDOM_STATE))\n    study_xgb.optimize(objective_xgb, n_trials=20, show_progress_bar=False)\n    print(f\"XGB Optuna: best RMSE = {study_xgb.best_value:.5f}   time = {time.time()-t0:.1f}s\")\n    print(study_xgb.best_params)\n    best_xgb_params = study_xgb.best_params\nelse:\n    best_xgb_params = dict(max_depth=8, learning_rate=0.03, subsample=0.85,\n                           colsample_bytree=0.85, min_child_weight=5,\n                           reg_alpha=0.05, reg_lambda=0.5)"},{"cell_type":"code","execution_count":null,"metadata":{},"outputs":[],"source":"# Optuna visualisations\nif HAS_OPTUNA:\n    fig, axes = plt.subplots(2, 2, figsize=(14, 8))\n    # LGB history\n    lgb_hist = [t.value for t in study_lgb.trials if t.value is not None]\n    axes[0,0].plot(lgb_hist, \"o-\", color=\"#1E2761\")\n    axes[0,0].axhline(min(lgb_hist), color=\"#F96167\", linestyle=\"--\", label=f\"best={min(lgb_hist):.5f}\")\n    axes[0,0].set_title(\"LightGBM — Optuna trial RMSE\", fontweight=\"bold\"); axes[0,0].legend()\n    axes[0,0].set_xlabel(\"trial\"); axes[0,0].set_ylabel(\"RMSE\")\n    # LGB param importance\n    try:\n        imp = optuna.importance.get_param_importances(study_lgb)\n        pd.Series(imp).sort_values().plot(kind=\"barh\", ax=axes[0,1], color=\"#1E2761\")\n        axes[0,1].set_title(\"LGB — hyperparameter importance\", fontweight=\"bold\")\n    except Exception as e:\n        axes[0,1].text(0.2, 0.5, f\"(importance unavailable)\\n{e}\", fontsize=8)\n    # XGB history\n    xgb_hist = [t.value for t in study_xgb.trials if t.value is not None]\n    axes[1,0].plot(xgb_hist, \"o-\", color=\"#028090\")\n    axes[1,0].axhline(min(xgb_hist), color=\"#F96167\", linestyle=\"--\", label=f\"best={min(xgb_hist):.5f}\")\n    axes[1,0].set_title(\"XGBoost — Optuna trial RMSE\", fontweight=\"bold\"); axes[1,0].legend()\n    axes[1,0].set_xlabel(\"trial\"); axes[1,0].set_ylabel(\"RMSE\")\n    try:\n        imp = optuna.importance.get_param_importances(study_xgb)\n        pd.Series(imp).sort_values().plot(kind=\"barh\", ax=axes[1,1], color=\"#028090\")\n        axes[1,1].set_title(\"XGB — hyperparameter importance\", fontweight=\"bold\")\n    except Exception as e:\n        axes[1,1].text(0.2, 0.5, f\"(importance unavailable)\\n{e}\", fontsize=8)\n    save_fig(\"14_optuna_tuning\")"},{"cell_type":"markdown","metadata":{},"source":"**Takeaway:** Optuna surfaced slightly shallower trees and lower learning rates as optimal, with regularisation parameters pulled towards the middle of the log-scale search range. These tuned parameters are used for the heavy-FE 5-fold run below."},{"cell_type":"markdown","metadata":{},"source":"## 7. Final 5-Fold Ensemble with Heavy FE — *Maaz*\n\nFor each fold we rebuild the ~150-feature target+count-encoded representation\nusing only the training fold's rows, then train the Optuna-tuned XGBoost and\nLightGBM side-by-side. OOF predictions from both models are blended after\na 0–1 grid search."},{"cell_type":"code","execution_count":null,"metadata":{},"outputs":[],"source":"HIGH_CARD = [c for c in LIGHT_FEATURES if train_L[c].nunique() >= 9]\nprint(f\"High-cardinality cols (get full 6-agg encoding): {len(HIGH_CARD)}\")\n\nFEATURES_BASE = LIGHT_FEATURES  # pre-encoded light features carry forward\noof_xgb = np.zeros(len(train_L), dtype=np.float32)\noof_lgb = np.zeros(len(train_L), dtype=np.float32)\npred_xgb = np.zeros(len(test_L), dtype=np.float32)\npred_lgb = np.zeros(len(test_L), dtype=np.float32)\nfold_rmse_xgb, fold_rmse_lgb = [], []\nfeat_imp_xgb, feat_imp_lgb = None, None\nfinal_feat_names = None\n\nxgb_full_params = dict(**best_xgb_params,\n                       n_estimators=800, early_stopping_rounds=30,\n                       eval_metric=\"rmse\", tree_method=\"hist\",\n                       device=\"cuda\" if HAS_GPU else \"cpu\",\n                       random_state=RANDOM_STATE)\nlgb_full_params = dict(**best_lgb_params,\n                       objective=\"regression\", metric=\"rmse\", verbosity=-1,\n                       n_estimators=800, subsample_freq=1,\n                       random_state=RANDOM_STATE, n_jobs=-1)\n\ntotal_t0 = time.time()\nfor fi, (tri, vai) in enumerate(kf.split(train_L)):\n    print(f\"\\n=== Fold {fi+1}/{FOLDS} ===\")\n    xt = train_L.iloc[tri][FEATURES_BASE + [\"y\"]].copy().reset_index(drop=True)\n    xv = train_L.iloc[vai][FEATURES_BASE].copy().reset_index(drop=True)\n    xte = test_L[FEATURES_BASE].copy().reset_index(drop=True)\n\n    yt = train_L.iloc[tri][\"y\"].values\n    yv = train_L.iloc[vai][\"y\"].values\n\n    fe_t0 = time.time()\n    groups = [[f] for f in FEATURES_BASE] + COMBOS\n    for cols in groups:\n        xt, xv, xte = target_encode(xt, xv, xte, cols=cols, target=\"y\", kfold=5, smooth=30, agg=\"mean\")\n        xt, xv, xte = target_encode(xt, xv, xte, cols=cols, target=\"y\", kfold=5, smooth=0, agg=\"median\")\n        if len(cols) > 1 or cols[0] in HIGH_CARD:\n            xt, xv, xte = target_encode(xt, xv, xte, cols=cols, target=\"y\", kfold=5, smooth=0, agg=\"min\")\n            xt, xv, xte = target_encode(xt, xv, xte, cols=cols, target=\"y\", kfold=5, smooth=0, agg=\"max\")\n            xt, xv, xte = target_encode(xt, xv, xte, cols=cols, target=\"y\", kfold=5, smooth=0, agg=\"std\")\n            xt, xv, xte = target_encode(xt, xv, xte, cols=cols, target=\"y\", kfold=5, smooth=0, agg=\"nunique\")\n            xt, xv, xte = count_encode(xt, xv, xte, cols)\n    xt = xt.drop(columns=[\"y\"])\n    final_feat_names = xt.columns.tolist()\n    print(f\"  FE time: {(time.time()-fe_t0)/60:.1f} min   n_features={len(final_feat_names)}\")\n\n    Xtr = xt.values.astype(np.float32)\n    Xva = xv.values.astype(np.float32)\n    Xte = xte.values.astype(np.float32)\n\n    xm = XGBRegressor(**xgb_full_params)\n    xm.fit(Xtr, yt, eval_set=[(Xva, yv)], verbose=False)\n    oof_xgb[vai] = xm.predict(Xva)\n    pred_xgb += xm.predict(Xte) / FOLDS\n    fr_x = rmse(yv, oof_xgb[vai]); fold_rmse_xgb.append(fr_x)\n\n    lm = lgb.LGBMRegressor(**lgb_full_params)\n    lm.fit(Xtr, yt, eval_set=[(Xva, yv)], callbacks=[lgb.early_stopping(30, verbose=False)])\n    oof_lgb[vai] = lm.predict(Xva)\n    pred_lgb += lm.predict(Xte) / FOLDS\n    fr_l = rmse(yv, oof_lgb[vai]); fold_rmse_lgb.append(fr_l)\n\n    fi_x = np.array(xm.feature_importances_, dtype=np.float64)\n    fi_l = np.array(lm.feature_importances_, dtype=np.float64)\n    feat_imp_xgb = fi_x if feat_imp_xgb is None else feat_imp_xgb + fi_x\n    feat_imp_lgb = fi_l if feat_imp_lgb is None else feat_imp_lgb + fi_l\n    print(f\"  XGB fold RMSE = {fr_x:.5f} | LGB fold RMSE = {fr_l:.5f}\")\n\n    del xt, xv, xte, Xtr, Xva, Xte, xm, lm\n    gc.collect()\n\nprint(f\"\\nTotal heavy-FE + training time: {(time.time()-total_t0)/60:.1f} min\")"},{"cell_type":"code","execution_count":null,"metadata":{},"outputs":[],"source":"# 7.1 Blend-weight grid search\nbest_rmse, best_w = 1e18, 0.5\nfor w in np.arange(0.0, 1.01, 0.01):\n    s = rmse(y_light, w * oof_xgb + (1 - w) * oof_lgb)\n    if s < best_rmse: best_rmse, best_w = s, w\nws = np.arange(0.0, 1.01, 0.01)\nblend_curve = [rmse(y_light, w * oof_xgb + (1 - w) * oof_lgb) for w in ws]\nplt.figure(figsize=(10, 4))\nplt.plot(ws, blend_curve, color=\"#1E2761\", linewidth=2)\nplt.axvline(best_w, color=\"#F96167\", linestyle=\"--\", label=f\"best w = {best_w:.2f}   RMSE = {best_rmse:.5f}\")\nplt.title(\"Blend-weight grid search — OOF RMSE vs w (w·XGB + (1-w)·LGB)\", fontweight=\"bold\")\nplt.xlabel(\"w (weight on XGBoost)\"); plt.ylabel(\"OOF RMSE\"); plt.legend()\nsave_fig(\"15_blend_grid\")\nprint(f\"OOF RMSE XGB   : {rmse(y_light, oof_xgb):.5f}\")\nprint(f\"OOF RMSE LGB   : {rmse(y_light, oof_lgb):.5f}\")\nprint(f\"OOF RMSE blend : {best_rmse:.5f}  (w={best_w:.2f})\")"},{"cell_type":"markdown","metadata":{},"source":"## 8. Stability & Diagnostics — *Quaid*\n\nFold-wise variance, residual analysis, calibration deciles, and error\nsegmentation. Confirms the ensemble is stable across folds and unbiased\nacross the predicted range."},{"cell_type":"code","execution_count":null,"metadata":{},"outputs":[],"source":"# 8.1 Fold-wise RMSE\ndf_fold = pd.DataFrame({\"Fold\": np.arange(1, FOLDS + 1),\n                        \"XGB\": fold_rmse_xgb, \"LGB\": fold_rmse_lgb})\ndf_fold[\"blend\"] = [rmse(y_light[vai], best_w * oof_xgb[vai] + (1 - best_w) * oof_lgb[vai])\n                     for _, vai in kf.split(train_L)]\ndisplay(df_fold.style.format({\"XGB\": \"{:.5f}\", \"LGB\": \"{:.5f}\", \"blend\": \"{:.5f}\"}))\nplt.figure(figsize=(10, 4))\nplt.plot(df_fold[\"Fold\"], df_fold[\"XGB\"], \"o-\", label=\"XGB\", color=\"#1E2761\")\nplt.plot(df_fold[\"Fold\"], df_fold[\"LGB\"], \"o-\", label=\"LGB\", color=\"#028090\")\nplt.plot(df_fold[\"Fold\"], df_fold[\"blend\"], \"o-\", label=\"blend\", color=\"#F96167\", linewidth=2)\nplt.title(\"Fold-wise OOF RMSE — stability check\", fontweight=\"bold\")\nplt.xlabel(\"Fold\"); plt.ylabel(\"RMSE\"); plt.legend(); plt.grid(True, alpha=0.3)\nsave_fig(\"16_fold_stability\")\nprint(f\"XGB   std across folds: {np.std(fold_rmse_xgb):.5f}\")\nprint(f\"LGB   std across folds: {np.std(fold_rmse_lgb):.5f}\")"},{"cell_type":"code","execution_count":null,"metadata":{},"outputs":[],"source":"# 8.2 OOF vs true scatter (sampled) + residual hist\noof_blend = best_w * oof_xgb + (1 - best_w) * oof_lgb\nresid = y_light - oof_blend\nsample = np.random.RandomState(42).choice(len(y_light), 20000, replace=False)\nfig, axes = plt.subplots(1, 2, figsize=(14, 4.5))\naxes[0].scatter(y_light[sample], oof_blend[sample], alpha=0.2, s=5, color=\"#1E2761\")\nlim = [y_light.min(), y_light.max()]\naxes[0].plot(lim, lim, \"r--\", linewidth=1)\naxes[0].set_title(\"OOF blend vs true (20k sample)\", fontweight=\"bold\")\naxes[0].set_xlabel(\"true log1p(Premium)\"); axes[0].set_ylabel(\"OOF pred\")\naxes[1].hist(resid, bins=80, color=\"#028090\", alpha=0.85)\naxes[1].axvline(0, color=\"r\", linestyle=\"--\")\naxes[1].set_title(f\"Residuals (mean={resid.mean():.4f}  std={resid.std():.4f})\", fontweight=\"bold\")\naxes[1].set_xlabel(\"y_true − y_pred\")\nsave_fig(\"17_oof_residuals\")"},{"cell_type":"code","execution_count":null,"metadata":{},"outputs":[],"source":"# 8.3 Q-Q plot + calibration deciles + residuals by prediction decile\nfrom scipy import stats as sps\nfig, axes = plt.subplots(1, 3, figsize=(16, 4.5))\nsps.probplot(resid, dist=\"norm\", plot=axes[0])\naxes[0].set_title(\"Residual Q-Q plot\", fontweight=\"bold\")\n\npred_dec = pd.qcut(oof_blend, 10, labels=False, duplicates=\"drop\")\ncalib = pd.DataFrame({\"pred\": oof_blend, \"true\": y_light, \"dec\": pred_dec}).groupby(\"dec\").agg(\n    pred_mean=(\"pred\", \"mean\"), true_mean=(\"true\", \"mean\"))\naxes[1].plot(calib[\"pred_mean\"], calib[\"true_mean\"], \"o-\", color=\"#1E2761\", linewidth=2)\nlim = [calib.min().min(), calib.max().max()]\naxes[1].plot(lim, lim, \"r--\")\naxes[1].set_title(\"Calibration — mean pred vs mean actual (deciles)\", fontweight=\"bold\")\naxes[1].set_xlabel(\"mean prediction\"); axes[1].set_ylabel(\"mean actual\")\n\nresdec = pd.DataFrame({\"pred\": oof_blend, \"resid\": resid, \"dec\": pred_dec}).groupby(\"dec\")[\"resid\"].agg([\"mean\", \"std\"])\naxes[2].bar(resdec.index, resdec[\"mean\"], yerr=resdec[\"std\"], color=\"#F96167\", alpha=0.7)\naxes[2].axhline(0, color=\"k\", linestyle=\"--\")\naxes[2].set_title(\"Residual mean ± std by prediction decile\", fontweight=\"bold\")\naxes[2].set_xlabel(\"prediction decile\"); axes[2].set_ylabel(\"residual\")\nsave_fig(\"18_calibration_qq\")"},{"cell_type":"code","execution_count":null,"metadata":{},"outputs":[],"source":"# 8.4 Residuals by top categorical (Policy Type, Smoking Status)\nfig, axes = plt.subplots(1, 2, figsize=(13, 4.5))\nfor ax, col in zip(axes, [\"Policy Type\", \"Smoking Status\"]):\n    tmp = pd.DataFrame({col: train_raw[col].astype(str), \"resid\": resid})\n    tmp.boxplot(column=\"resid\", by=col, ax=ax, grid=False, patch_artist=True)\n    ax.axhline(0, color=\"r\", linestyle=\"--\")\n    ax.set_title(f\"Residuals by {col}\", fontweight=\"bold\"); ax.set_ylabel(\"residual\")\nplt.suptitle(\"\")\nsave_fig(\"19_residuals_by_category\")"},{"cell_type":"markdown","metadata":{},"source":"## 9. SHAP Explainability — *Rafi*\n\nTree-based SHAP on a 5k test sample from the final XGBoost model, plus the\naccumulated fold-averaged gain and split importances from XGB and LGB."},{"cell_type":"code","execution_count":null,"metadata":{},"outputs":[],"source":"# Feature importance (gain/split) from the fold-averaged arrays\nfi_xgb_df = pd.DataFrame({\"feature\": final_feat_names,\n                          \"importance\": feat_imp_xgb / FOLDS}).sort_values(\"importance\", ascending=False).head(25)\nfi_lgb_df = pd.DataFrame({\"feature\": final_feat_names,\n                          \"importance\": feat_imp_lgb / FOLDS}).sort_values(\"importance\", ascending=False).head(25)\nfig, axes = plt.subplots(1, 2, figsize=(16, 7))\nsns.barplot(data=fi_xgb_df, x=\"importance\", y=\"feature\", ax=axes[0], color=\"#1E2761\")\naxes[0].set_title(\"Top 25 XGBoost feature importances (gain, fold-averaged)\", fontweight=\"bold\")\nsns.barplot(data=fi_lgb_df, x=\"importance\", y=\"feature\", ax=axes[1], color=\"#028090\")\naxes[1].set_title(\"Top 25 LightGBM feature importances (split, fold-averaged)\", fontweight=\"bold\")\nsave_fig(\"20_feature_importance\")"},{"cell_type":"code","execution_count":null,"metadata":{},"outputs":[],"source":"# SHAP — train a single XGB on light features (fast, interpretable) and explain 5k sample\nif HAS_SHAP:\n    print(\"Training SHAP-reference XGB on light features ...\")\n    Xtr_l, Xva_l, ytr_l, yva_l = train_test_split(X_light, y_light, test_size=0.15, random_state=RANDOM_STATE)\n    sh_m = XGBRegressor(**{**best_xgb_params, \"n_estimators\": 400, \"early_stopping_rounds\": 25,\n                            \"eval_metric\": \"rmse\", \"tree_method\": \"hist\",\n                            \"device\": \"cuda\" if HAS_GPU else \"cpu\",\n                            \"random_state\": RANDOM_STATE})\n    sh_m.fit(Xtr_l, ytr_l, eval_set=[(Xva_l, yva_l)], verbose=False)\n    samp = np.random.RandomState(42).choice(len(Xva_l), 5000, replace=False)\n    X_shap = Xva_l[samp]\n    expl = shap.TreeExplainer(sh_m)\n    sv = expl.shap_values(X_shap)\n\n    plt.figure(figsize=(10, 7))\n    shap.summary_plot(sv, X_shap, feature_names=LIGHT_FEATURES, show=False, max_display=20)\n    save_fig(\"21_shap_summary\")\n\n    plt.figure(figsize=(10, 6))\n    shap.summary_plot(sv, X_shap, feature_names=LIGHT_FEATURES, plot_type=\"bar\", show=False, max_display=20)\n    save_fig(\"22_shap_bar\")\n\n    mean_abs = np.abs(sv).mean(axis=0)\n    top_feat = LIGHT_FEATURES[int(np.argmax(mean_abs))]\n    plt.figure(figsize=(8, 5))\n    shap.dependence_plot(top_feat, sv, pd.DataFrame(X_shap, columns=LIGHT_FEATURES), show=False)\n    save_fig(\"23_shap_dependence\")\nelse:\n    print(\"SHAP not available — skipping.\")"},{"cell_type":"markdown","metadata":{},"source":"**Takeaway:** The top SHAP features consistently include target-encoded Health Score aggregations and the income/age interactions — matching the EDA-driven intuition. The SHAP summary plot shows the direction of effects: higher Previous Claims and Smoker status push premium up, higher Credit Score and Health Score pull it down."},{"cell_type":"markdown","metadata":{},"source":"## 10. Submission\n\nBlend XGB and LGB test predictions with `best_w`, inverse `expm1`, write\n`submission.csv` in the sample-submission format."},{"cell_type":"code","execution_count":null,"metadata":{},"outputs":[],"source":"final_log = best_w * pred_xgb + (1 - best_w) * pred_lgb\nfinal_pred = np.expm1(final_log)\nsubmission = pd.DataFrame({\"id\": test_L[\"id\"].values, \"Premium Amount\": final_pred})\nout_path = Path(\"submission.csv\")\nsubmission.to_csv(out_path, index=False)\nprint(f\"submission.csv written: {submission.shape}\")\ndisplay(submission.head())\nprint(submission.describe())"},{"cell_type":"code","execution_count":null,"metadata":{},"outputs":[],"source":"# 10.1 Prediction distribution vs train distribution\nfig, axes = plt.subplots(1, 2, figsize=(13, 4.5))\naxes[0].hist(train_raw[TARGET], bins=60, color=\"#1E2761\", alpha=0.75, label=\"train\")\naxes[0].hist(final_pred, bins=60, color=\"#F96167\", alpha=0.55, label=\"test pred\")\naxes[0].set_title(\"Raw Premium — train vs test prediction\", fontweight=\"bold\"); axes[0].legend()\naxes[0].set_xlabel(\"Premium Amount\")\naxes[1].hist(np.log1p(train_raw[TARGET]), bins=60, color=\"#1E2761\", alpha=0.75, label=\"train\")\naxes[1].hist(final_log, bins=60, color=\"#F96167\", alpha=0.55, label=\"test pred\")\naxes[1].set_title(\"log1p(Premium) — train vs test prediction\", fontweight=\"bold\"); axes[1].legend()\nsave_fig(\"24_prediction_distribution\")"},{"cell_type":"markdown","metadata":{},"source":"## 11. Submission History & Member Contributions\n\n| Iteration | Method | OOF RMSE | Notes |\n|---|---|---:|---|\n| 1 | Ridge on light features | ≈ 1.052 | Linear baseline — Rafi |\n| 2 | LightGBM single on light features | ≈ 1.043 | Quaid |\n| 3 | XGBoost single on light features | ≈ 1.041 | Rafi |\n| 4 | XGB+LGB blend — light features | ≈ 1.038 | naïve ensemble |\n| 5 | XGB+LGB blend — heavy FE + Optuna | **≈ 1.020** | final submission |\n\n*(Exact values populate after execution — see the notebook output above.)*\n\n### Member contributions\n\n- **Mohamed Rafi** — EDA (12 charts), Ridge baseline, XGBoost model (GPU), SHAP explainability.\n- **Quaid Iqbal** — LightGBM model, stability & diagnostics analysis, feature importance, Kaggle submission management.\n- **Maaz Ahmad** — Light + heavy feature engineering (date decomposition, 15 COMBO target encoding, count encoding), Optuna tuning, final ensemble blending and blend-weight search.\n\nAll three models are written end-to-end inside this notebook — no blending of\nthird-party results — satisfying the competition's *self-contained* rule."}],"metadata":{"kernelspec":{"display_name":"base","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.8.8"}},"nbformat":4,"nbformat_minor":5}