Python 64.6%
TeX 23.3%
HTML 11.6%
Makefile 0.5%
1"""2ml_models.py3============4Machine learning models for hedonic pricing analysis of Airbnb listings.56Provides training, evaluation, cross-validation, feature importance extraction,7and hyperparameter tuning for a battery of regression models including8ElasticNet, Random Forest, Gradient Boosting, XGBoost, LightGBM, CatBoost,9and Generalized Additive Models (GAM).1011Author: Simon-Pierre Boucher12Date: 2026-05-2013"""1415from __future__ import annotations1617import logging18import warnings19from typing import Any, Dict, List, Optional, Tuple, Union2021import numpy as np22import pandas as pd23from sklearn.ensemble import GradientBoostingRegressor, RandomForestRegressor24from sklearn.linear_model import ElasticNetCV25from sklearn.metrics import (26 mean_absolute_error,27 mean_absolute_percentage_error,28 mean_squared_error,29 r2_score,30)31from sklearn.model_selection import (32 KFold,33 RandomizedSearchCV,34 cross_val_score,35 train_test_split,36)37from sklearn.pipeline import Pipeline38from sklearn.preprocessing import PolynomialFeatures, StandardScaler3940logger = logging.getLogger(__name__)4142# ---------------------------------------------------------------------------43# Data Preparation44# ---------------------------------------------------------------------------454647def prepare_ml_data(48 df: pd.DataFrame,49 y_col: str,50 x_cols: List[str],51 test_size: float = 0.2,52 random_state: int = 42,53) -> Tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray, List[str]]:54 """Split data into train/test sets with NaN handling.5556 Parameters57 ----------58 df : pd.DataFrame59 Source dataframe.60 y_col : str61 Name of the target (price) column.62 x_cols : list[str]63 Feature column names.64 test_size : float65 Fraction held out for testing.66 random_state : int67 Seed for reproducibility.6869 Returns70 -------71 X_train, X_test, y_train, y_test : np.ndarray72 Split arrays.73 feature_names : list[str]74 Column names retained after NaN removal.75 """76 subset = df[x_cols + [y_col]].copy()7778 n_before = len(subset)79 subset = subset.dropna()80 n_after = len(subset)81 if n_before != n_after:82 logger.info(83 "Dropped %d rows containing NaN (%d -> %d).",84 n_before - n_after,85 n_before,86 n_after,87 )8889 X = subset[x_cols].values.astype(np.float64)90 y = subset[y_col].values.astype(np.float64)91 feature_names = list(x_cols)9293 X_train, X_test, y_train, y_test = train_test_split(94 X, y, test_size=test_size, random_state=random_state95 )9697 logger.info(98 "Train: %d samples, Test: %d samples, Features: %d",99 X_train.shape[0],100 X_test.shape[0],101 X_train.shape[1],102 )103 return X_train, X_test, y_train, y_test, feature_names104105106# ---------------------------------------------------------------------------107# Individual Model Trainers108# ---------------------------------------------------------------------------109110111def train_elastic_net(112 X_train: np.ndarray,113 y_train: np.ndarray,114 alphas: Optional[List[float]] = None,115 l1_ratios: Optional[List[float]] = None,116 cv: int = 5,117 random_state: int = 42,118) -> ElasticNetCV:119 """Train an ElasticNet with built-in cross-validated alpha selection.120121 Parameters122 ----------123 X_train, y_train : np.ndarray124 Training data.125 alphas : list[float] | None126 Candidate regularisation strengths. Defaults to a log-spaced grid.127 l1_ratios : list[float] | None128 L1/L2 mixing ratios to try. Defaults to ``[0.1, 0.5, 0.7, 0.9, 0.95, 1.0]``.129 cv : int130 Number of CV folds.131 random_state : int132 Seed for reproducibility.133134 Returns135 -------136 ElasticNetCV137 Fitted model.138 """139 if alphas is None:140 alphas = np.logspace(-4, 2, 50).tolist()141 if l1_ratios is None:142 l1_ratios = [0.1, 0.5, 0.7, 0.9, 0.95, 1.0]143144 model = ElasticNetCV(145 alphas=alphas,146 l1_ratio=l1_ratios,147 cv=cv,148 random_state=random_state,149 max_iter=10_000,150 n_jobs=-1,151 )152 model.fit(X_train, y_train)153 logger.info(154 "ElasticNet alpha=%.6f l1_ratio=%.2f",155 model.alpha_,156 model.l1_ratio_,157 )158 return model159160161def train_random_forest(162 X_train: np.ndarray,163 y_train: np.ndarray,164 n_estimators: int = 500,165 max_depth: Optional[int] = None,166 min_samples_leaf: int = 5,167 random_state: int = 42,168) -> RandomForestRegressor:169 """Train a Random Forest regressor.170171 Parameters172 ----------173 X_train, y_train : np.ndarray174 Training data.175 n_estimators : int176 Number of trees.177 max_depth : int | None178 Maximum tree depth (``None`` = unlimited).179 min_samples_leaf : int180 Minimum samples per leaf node.181 random_state : int182 Seed for reproducibility.183184 Returns185 -------186 RandomForestRegressor187 Fitted model.188 """189 model = RandomForestRegressor(190 n_estimators=n_estimators,191 max_depth=max_depth,192 min_samples_leaf=min_samples_leaf,193 random_state=random_state,194 n_jobs=-1,195 )196 model.fit(X_train, y_train)197 logger.info("RandomForest n_estimators=%d fitted.", n_estimators)198 return model199200201def train_gradient_boosting(202 X_train: np.ndarray,203 y_train: np.ndarray,204 n_estimators: int = 500,205 learning_rate: float = 0.05,206 max_depth: int = 5,207 subsample: float = 0.8,208 random_state: int = 42,209) -> GradientBoostingRegressor:210 """Train a scikit-learn GradientBoostingRegressor.211212 Parameters213 ----------214 X_train, y_train : np.ndarray215 Training data.216 n_estimators : int217 Boosting rounds.218 learning_rate : float219 Shrinkage factor.220 max_depth : int221 Maximum tree depth per stage.222 subsample : float223 Row subsampling fraction.224 random_state : int225 Seed for reproducibility.226227 Returns228 -------229 GradientBoostingRegressor230 Fitted model.231 """232 model = GradientBoostingRegressor(233 n_estimators=n_estimators,234 learning_rate=learning_rate,235 max_depth=max_depth,236 subsample=subsample,237 random_state=random_state,238 )239 model.fit(X_train, y_train)240 logger.info("GradientBoosting fitted (n_estimators=%d).", n_estimators)241 return model242243244def train_xgboost(245 X_train: np.ndarray,246 y_train: np.ndarray,247 n_estimators: int = 500,248 learning_rate: float = 0.05,249 max_depth: int = 6,250 subsample: float = 0.8,251 colsample_bytree: float = 0.8,252 random_state: int = 42,253) -> Any:254 """Train an XGBRegressor. Falls back to sklearn GradientBoostingRegressor255 if ``xgboost`` is not installed.256257 Returns258 -------259 model260 Fitted regressor (XGBRegressor or GradientBoostingRegressor).261 """262 try:263 from xgboost import XGBRegressor264265 model = XGBRegressor(266 n_estimators=n_estimators,267 learning_rate=learning_rate,268 max_depth=max_depth,269 subsample=subsample,270 colsample_bytree=colsample_bytree,271 random_state=random_state,272 n_jobs=-1,273 verbosity=0,274 )275 model.fit(X_train, y_train)276 logger.info("XGBoost fitted (n_estimators=%d).", n_estimators)277 return model278 except ImportError:279 logger.warning("xgboost not installed — falling back to sklearn GBR.")280 return train_gradient_boosting(281 X_train,282 y_train,283 n_estimators=n_estimators,284 learning_rate=learning_rate,285 max_depth=max_depth,286 subsample=subsample,287 random_state=random_state,288 )289290291def train_lightgbm(292 X_train: np.ndarray,293 y_train: np.ndarray,294 n_estimators: int = 500,295 learning_rate: float = 0.05,296 num_leaves: int = 31,297 subsample: float = 0.8,298 colsample_bytree: float = 0.8,299 random_state: int = 42,300) -> Any:301 """Train an LGBMRegressor. Falls back to sklearn GBR if ``lightgbm`` is302 not installed.303304 Returns305 -------306 model307 Fitted regressor.308 """309 try:310 from lightgbm import LGBMRegressor311312 model = LGBMRegressor(313 n_estimators=n_estimators,314 learning_rate=learning_rate,315 num_leaves=num_leaves,316 subsample=subsample,317 colsample_bytree=colsample_bytree,318 random_state=random_state,319 n_jobs=-1,320 verbose=-1,321 )322 model.fit(X_train, y_train)323 logger.info("LightGBM fitted (n_estimators=%d).", n_estimators)324 return model325 except ImportError:326 logger.warning("lightgbm not installed — falling back to sklearn GBR.")327 return train_gradient_boosting(328 X_train,329 y_train,330 n_estimators=n_estimators,331 learning_rate=learning_rate,332 max_depth=5,333 subsample=subsample,334 random_state=random_state,335 )336337338def train_catboost(339 X_train: np.ndarray,340 y_train: np.ndarray,341 iterations: int = 500,342 learning_rate: float = 0.05,343 depth: int = 6,344 random_state: int = 42,345) -> Any:346 """Train a CatBoostRegressor. Falls back to sklearn GBR if ``catboost``347 is not installed.348349 Returns350 -------351 model352 Fitted regressor.353 """354 try:355 from catboost import CatBoostRegressor356357 model = CatBoostRegressor(358 iterations=iterations,359 learning_rate=learning_rate,360 depth=depth,361 random_state=random_state,362 verbose=0,363 )364 model.fit(X_train, y_train)365 logger.info("CatBoost fitted (iterations=%d).", iterations)366 return model367 except ImportError:368 logger.warning("catboost not installed — falling back to sklearn GBR.")369 return train_gradient_boosting(370 X_train,371 y_train,372 n_estimators=iterations,373 learning_rate=learning_rate,374 max_depth=depth,375 random_state=random_state,376 )377378379def train_gam(380 X_train: np.ndarray,381 y_train: np.ndarray,382 n_splines: int = 20,383 degree: int = 2,384) -> Any:385 """Train a Generalized Additive Model.386387 Attempts to use ``pygam.LinearGAM``; if unavailable, builds a polynomial-388 feature pipeline as a lightweight approximation.389390 Returns391 -------392 model393 Fitted GAM or polynomial pipeline.394 """395 try:396 from pygam import LinearGAM, s397398 # Build one spline term per feature.399 terms = s(0, n_splines=n_splines)400 for j in range(1, X_train.shape[1]):401 terms += s(j, n_splines=n_splines)402403 model = LinearGAM(terms)404 model.gridsearch(X_train, y_train, progress=False)405 logger.info("pygam LinearGAM fitted (n_splines=%d).", n_splines)406 return model407 except ImportError:408 logger.warning(409 "pygam not installed — approximating GAM with polynomial features (degree=%d).",410 degree,411 )412 pipe = Pipeline(413 [414 ("poly", PolynomialFeatures(degree=degree, interaction_only=True, include_bias=False)),415 ("scaler", StandardScaler()),416 (417 "enet",418 ElasticNetCV(419 l1_ratio=[0.5, 0.7, 0.9, 0.95, 1.0],420 cv=5,421 max_iter=10_000,422 n_jobs=-1,423 ),424 ),425 ]426 )427 pipe.fit(X_train, y_train)428 logger.info("Polynomial GAM approximation fitted (degree=%d).", degree)429 return pipe430431432# ---------------------------------------------------------------------------433# Train All Models434# ---------------------------------------------------------------------------435436437def train_all_models(438 X_train: np.ndarray,439 y_train: np.ndarray,440) -> Dict[str, Any]:441 """Train every available model and return a name -> fitted model dict.442443 Models that fail to train (e.g. convergence issues) are logged and444 skipped rather than raising.445446 Returns447 -------448 dict[str, model]449 Mapping of model name to fitted estimator.450 """451 trainers: Dict[str, Any] = {452 "ElasticNet": train_elastic_net,453 "RandomForest": train_random_forest,454 "GradientBoosting": train_gradient_boosting,455 "XGBoost": train_xgboost,456 "LightGBM": train_lightgbm,457 "CatBoost": train_catboost,458 "GAM": train_gam,459 }460461 models: Dict[str, Any] = {}462 for name, trainer in trainers.items():463 try:464 logger.info("Training %s ...", name)465 models[name] = trainer(X_train, y_train)466 except Exception:467 logger.exception("Failed to train %s — skipping.", name)468 return models469470471# ---------------------------------------------------------------------------472# Evaluation473# ---------------------------------------------------------------------------474475476def evaluate_model(477 model: Any,478 X_test: np.ndarray,479 y_test: np.ndarray,480) -> Dict[str, float]:481 """Compute regression metrics for a fitted model.482483 Returns484 -------485 dict486 Keys: ``R2``, ``RMSE``, ``MAE``, ``MAPE``.487 """488 y_pred = model.predict(X_test)489 r2 = r2_score(y_test, y_pred)490 rmse = float(np.sqrt(mean_squared_error(y_test, y_pred)))491 mae = mean_absolute_error(y_test, y_pred)492493 # Guard against zero prices when computing MAPE.494 with warnings.catch_warnings():495 warnings.simplefilter("ignore")496 mask = y_test != 0497 if mask.sum() > 0:498 mape = float(np.mean(np.abs((y_test[mask] - y_pred[mask]) / y_test[mask]))) * 100499 else:500 mape = np.nan501502 metrics = {"R2": r2, "RMSE": rmse, "MAE": mae, "MAPE": mape}503 logger.info(504 "Evaluation R2=%.4f RMSE=%.2f MAE=%.2f MAPE=%.2f%%",505 r2,506 rmse,507 mae,508 mape,509 )510 return metrics511512513def compare_all_models(514 models_dict: Dict[str, Any],515 X_test: np.ndarray,516 y_test: np.ndarray,517) -> pd.DataFrame:518 """Evaluate every model and return a comparison DataFrame sorted by R2.519520 Parameters521 ----------522 models_dict : dict[str, model]523 Output of :func:`train_all_models`.524 X_test, y_test : np.ndarray525 Held-out test data.526527 Returns528 -------529 pd.DataFrame530 Columns: ``Model``, ``R2``, ``RMSE``, ``MAE``, ``MAPE``.531 """532 rows: List[Dict[str, Any]] = []533 for name, model in models_dict.items():534 try:535 metrics = evaluate_model(model, X_test, y_test)536 rows.append({"Model": name, **metrics})537 except Exception:538 logger.exception("Evaluation failed for %s.", name)539540 comparison = pd.DataFrame(rows).sort_values("R2", ascending=False).reset_index(drop=True)541 logger.info("\n%s", comparison.to_string(index=False))542 return comparison543544545# ---------------------------------------------------------------------------546# Cross-Validation547# ---------------------------------------------------------------------------548549550def cross_validate_model(551 model: Any,552 X: np.ndarray,553 y: np.ndarray,554 k: int = 5,555 random_state: int = 42,556) -> Dict[str, float]:557 """K-fold cross-validation returning mean and std of R2, RMSE, and MAE.558559 Parameters560 ----------561 model : estimator562 An **unfitted** or fitted sklearn-compatible estimator. A clone is563 created internally, so the caller's model is not mutated.564 X, y : np.ndarray565 Full dataset (train + test).566 k : int567 Number of folds.568 random_state : int569 Seed for KFold split.570571 Returns572 -------573 dict574 Keys: ``R2_mean``, ``R2_std``, ``RMSE_mean``, ``RMSE_std``,575 ``MAE_mean``, ``MAE_std``.576 """577 from sklearn.base import clone578579 kf = KFold(n_splits=k, shuffle=True, random_state=random_state)580 r2_scores: List[float] = []581 rmse_scores: List[float] = []582 mae_scores: List[float] = []583584 for train_idx, val_idx in kf.split(X):585 X_tr, X_val = X[train_idx], X[val_idx]586 y_tr, y_val = y[train_idx], y[val_idx]587588 m = clone(model)589 m.fit(X_tr, y_tr)590 y_pred = m.predict(X_val)591592 r2_scores.append(r2_score(y_val, y_pred))593 rmse_scores.append(float(np.sqrt(mean_squared_error(y_val, y_pred))))594 mae_scores.append(mean_absolute_error(y_val, y_pred))595596 results = {597 "R2_mean": float(np.mean(r2_scores)),598 "R2_std": float(np.std(r2_scores)),599 "RMSE_mean": float(np.mean(rmse_scores)),600 "RMSE_std": float(np.std(rmse_scores)),601 "MAE_mean": float(np.mean(mae_scores)),602 "MAE_std": float(np.std(mae_scores)),603 }604 logger.info(605 "CV(%d) R2=%.4f +/- %.4f RMSE=%.2f +/- %.2f",606 k,607 results["R2_mean"],608 results["R2_std"],609 results["RMSE_mean"],610 results["RMSE_std"],611 )612 return results613614615# ---------------------------------------------------------------------------616# Feature Importance617# ---------------------------------------------------------------------------618619620def get_feature_importance(621 model: Any,622 feature_names: List[str],623) -> pd.DataFrame:624 """Extract feature importances from a fitted model.625626 Supports tree-based models (``feature_importances_``), linear models627 (``coef_``), and ``pygam`` (partial dependence variance).628629 Parameters630 ----------631 model : estimator632 Fitted model.633 feature_names : list[str]634 Feature column names.635636 Returns637 -------638 pd.DataFrame639 Columns: ``Feature``, ``Importance``, sorted descending.640 """641 importances: Optional[np.ndarray] = None642643 # Tree-based models (RF, GBR, XGB, LGBM, CatBoost).644 if hasattr(model, "feature_importances_"):645 importances = np.asarray(model.feature_importances_)646647 # Linear models (ElasticNet, Ridge, Lasso).648 elif hasattr(model, "coef_"):649 importances = np.abs(np.asarray(model.coef_))650651 # pygam — approximate importance via partial dependence variance.652 elif hasattr(model, "partial_dependence"):653 try:654 variances = []655 for j in range(len(feature_names)):656 XX = model.generate_X_grid(term=j)657 pdep, _ = model.partial_dependence(term=j, X=XX, width=0.95)658 variances.append(float(np.var(pdep)))659 importances = np.array(variances)660 except Exception:661 logger.warning("Could not extract importance from pygam model.")662663 # Pipeline fallback — try last step.664 elif hasattr(model, "named_steps"):665 last_step = list(model.named_steps.values())[-1]666 if hasattr(last_step, "coef_"):667 # Polynomial pipeline: importances don't map 1-to-1 to original668 # features. Aggregate absolute coefficients.669 poly = model.named_steps.get("poly")670 if poly is not None:671 coefs = np.abs(last_step.coef_)672 powers = poly.powers_673 agg = np.zeros(len(feature_names))674 for idx, power_vec in enumerate(powers):675 for feat_idx in np.nonzero(power_vec)[0]:676 if feat_idx < len(feature_names):677 agg[feat_idx] += coefs[idx] if idx < len(coefs) else 0.0678 importances = agg679 else:680 importances = np.abs(np.asarray(last_step.coef_))681682 if importances is None:683 logger.warning(684 "Cannot extract importance from model type %s.",685 type(model).__name__,686 )687 return pd.DataFrame(columns=["Feature", "Importance"])688689 # Truncate / pad to match feature_names length.690 if len(importances) > len(feature_names):691 importances = importances[: len(feature_names)]692 elif len(importances) < len(feature_names):693 importances = np.pad(importances, (0, len(feature_names) - len(importances)))694695 df_imp = (696 pd.DataFrame({"Feature": feature_names, "Importance": importances})697 .sort_values("Importance", ascending=False)698 .reset_index(drop=True)699 )700 return df_imp701702703# ---------------------------------------------------------------------------704# Hyperparameter Tuning705# ---------------------------------------------------------------------------706707708def tune_best_model(709 X_train: np.ndarray,710 y_train: np.ndarray,711 model_type: str = "xgboost",712 n_iter: int = 50,713 cv: int = 5,714 random_state: int = 42,715 scoring: str = "neg_root_mean_squared_error",716) -> Tuple[Any, Dict[str, Any]]:717 """Randomised hyperparameter search for a given model type.718719 Parameters720 ----------721 X_train, y_train : np.ndarray722 Training data.723 model_type : str724 One of ``'xgboost'``, ``'lightgbm'``, ``'catboost'``,725 ``'random_forest'``, ``'gradient_boosting'``.726 n_iter : int727 Number of parameter settings sampled.728 cv : int729 Cross-validation folds.730 random_state : int731 Seed for reproducibility.732 scoring : str733 Metric for ``RandomizedSearchCV``.734735 Returns736 -------737 best_model738 Refitted model with the best parameter set.739 best_params : dict740 Best parameters found.741 """742743 # --- estimator + param distributions -----------------------------------744 if model_type == "xgboost":745 try:746 from xgboost import XGBRegressor747748 estimator = XGBRegressor(random_state=random_state, n_jobs=-1, verbosity=0)749 except ImportError:750 logger.warning("xgboost unavailable — tuning sklearn GBR instead.")751 estimator = GradientBoostingRegressor(random_state=random_state)752 model_type = "gradient_boosting"753 elif model_type == "lightgbm":754 try:755 from lightgbm import LGBMRegressor756757 estimator = LGBMRegressor(random_state=random_state, n_jobs=-1, verbose=-1)758 except ImportError:759 logger.warning("lightgbm unavailable — tuning sklearn GBR instead.")760 estimator = GradientBoostingRegressor(random_state=random_state)761 model_type = "gradient_boosting"762 elif model_type == "catboost":763 try:764 from catboost import CatBoostRegressor765766 estimator = CatBoostRegressor(random_state=random_state, verbose=0)767 except ImportError:768 logger.warning("catboost unavailable — tuning sklearn GBR instead.")769 estimator = GradientBoostingRegressor(random_state=random_state)770 model_type = "gradient_boosting"771 elif model_type == "random_forest":772 estimator = RandomForestRegressor(random_state=random_state, n_jobs=-1)773 else:774 estimator = GradientBoostingRegressor(random_state=random_state)775 model_type = "gradient_boosting"776777 param_grids: Dict[str, Dict[str, Any]] = {778 "xgboost": {779 "n_estimators": [200, 500, 800, 1000],780 "learning_rate": [0.01, 0.03, 0.05, 0.1],781 "max_depth": [3, 4, 5, 6, 8],782 "subsample": [0.6, 0.7, 0.8, 0.9, 1.0],783 "colsample_bytree": [0.6, 0.7, 0.8, 0.9, 1.0],784 "min_child_weight": [1, 3, 5, 7],785 "reg_alpha": [0, 0.01, 0.1, 1],786 "reg_lambda": [0.5, 1, 2, 5],787 },788 "lightgbm": {789 "n_estimators": [200, 500, 800, 1000],790 "learning_rate": [0.01, 0.03, 0.05, 0.1],791 "num_leaves": [15, 31, 63, 127],792 "subsample": [0.6, 0.7, 0.8, 0.9, 1.0],793 "colsample_bytree": [0.6, 0.7, 0.8, 0.9, 1.0],794 "min_child_samples": [5, 10, 20, 50],795 "reg_alpha": [0, 0.01, 0.1, 1],796 "reg_lambda": [0.5, 1, 2, 5],797 },798 "catboost": {799 "iterations": [200, 500, 800, 1000],800 "learning_rate": [0.01, 0.03, 0.05, 0.1],801 "depth": [4, 6, 8, 10],802 "l2_leaf_reg": [1, 3, 5, 7, 9],803 "bagging_temperature": [0, 0.5, 1],804 },805 "random_forest": {806 "n_estimators": [200, 500, 800, 1000],807 "max_depth": [None, 10, 20, 30, 50],808 "min_samples_leaf": [1, 2, 5, 10],809 "min_samples_split": [2, 5, 10],810 "max_features": ["sqrt", "log2", 0.5, 0.8, None],811 },812 "gradient_boosting": {813 "n_estimators": [200, 500, 800, 1000],814 "learning_rate": [0.01, 0.03, 0.05, 0.1],815 "max_depth": [3, 4, 5, 6, 8],816 "subsample": [0.6, 0.7, 0.8, 0.9, 1.0],817 "min_samples_leaf": [1, 2, 5, 10],818 },819 }820821 param_dist = param_grids.get(model_type, param_grids["gradient_boosting"])822823 search = RandomizedSearchCV(824 estimator,825 param_distributions=param_dist,826 n_iter=n_iter,827 scoring=scoring,828 cv=cv,829 random_state=random_state,830 n_jobs=-1,831 verbose=0,832 )833 search.fit(X_train, y_train)834835 logger.info(836 "Best %s params: %s (score=%.4f)",837 model_type,838 search.best_params_,839 search.best_score_,840 )841 return search.best_estimator_, search.best_params_842