ramseshk
11182b47f8
Calibration overhaul: - Logistic regression mode for synthetic/bootstrap data (prevents LightGBM overfit) - 3-layer calibration stack: raw LR → Platt scaling → isotonic regression - Extreme probability smoothing: blend toward 0.5 when raw>0.95 or raw<0.05 - Platt preferred over isotonic (isotonic produces step functions with few points) - Continuous precipitation probability in bootstrap (beta distribution, not just 0/100) - Realistic NWP forecast errors: temp σ=2.0°C, rain calibration bias, diurnal-aware noise - Outlier injection: 10% of days have 2-3x larger errors (typhoon/low-pressure days) - LR model + StandardScaler saved as _lr.pkl alongside .lgb marker Results: - temp_gt_30c: AUC=0.987, Brier=0.049, predictions vary 20-85% per day - rain_gt_0mm: AUC=0.979, Brier=0.042, predictions vary 15-85% per day - temp_gt_35c: AUC=0.713 (realistic — extreme heat is hard to predict)
443 lines
16 KiB
Python
443 lines
16 KiB
Python
"""LightGBM probability models for HK weather prediction targets.
|
|
|
|
Each model uses a 3-layer calibration stack:
|
|
Layer 1: LightGBM binary classifier → raw log-odds
|
|
Layer 2: Platt scaling (logistic regression on validation logits)
|
|
Layer 3: Isotonic regression fallback (non-linear calibration)
|
|
|
|
Calibration parameters are saved/loaded with each model.
|
|
"""
|
|
|
|
import os
|
|
import json
|
|
import pickle
|
|
from pathlib import Path
|
|
from typing import Dict, Optional, Tuple, List
|
|
|
|
import numpy as np
|
|
import pandas as pd
|
|
|
|
try:
|
|
import lightgbm as lgb
|
|
except ImportError:
|
|
lgb = None
|
|
|
|
try:
|
|
from sklearn.isotonic import IsotonicRegression
|
|
from sklearn.linear_model import LogisticRegression
|
|
except ImportError:
|
|
IsotonicRegression = None
|
|
LogisticRegression = None
|
|
|
|
from config import DATA_DIR, PROJECT_ROOT
|
|
|
|
|
|
MODEL_DIR = Path(DATA_DIR) / "models"
|
|
|
|
|
|
TARGET_DEFINITIONS = {
|
|
"rain_gt_0mm_24h": {
|
|
"variable": "precipitation_sum", "threshold": 0.0, "op": "gt",
|
|
"description": "Precipitation > 0mm at t+24h",
|
|
},
|
|
"rain_gt_5mm_24h": {
|
|
"variable": "precipitation_sum", "threshold": 5.0, "op": "gt",
|
|
"description": "Precipitation > 5mm at t+24h",
|
|
},
|
|
"rain_gt_10mm_24h": {
|
|
"variable": "precipitation_sum", "threshold": 10.0, "op": "gt",
|
|
"description": "Precipitation > 10mm at t+24h",
|
|
},
|
|
"temp_gt_30c_24h": {
|
|
"variable": "temperature_2m_max", "threshold": 30.0, "op": "gt",
|
|
"description": "Tmax > 30°C at t+24h",
|
|
},
|
|
"temp_gt_33c_24h": {
|
|
"variable": "temperature_2m_max", "threshold": 33.0, "op": "gt",
|
|
"description": "Tmax > 33°C at t+24h",
|
|
},
|
|
"temp_gt_35c_24h": {
|
|
"variable": "temperature_2m_max", "threshold": 35.0, "op": "gt",
|
|
"description": "Tmax > 35°C at t+24h",
|
|
},
|
|
"wind_gt_30kmh_24h": {
|
|
"variable": "wind_speed_10m_max", "threshold": 30.0, "op": "gt",
|
|
"description": "Wind gust > 30 km/h at t+24h",
|
|
},
|
|
}
|
|
|
|
LGBM_PARAMS = {
|
|
"objective": "binary",
|
|
"metric": "binary_logloss",
|
|
"boosting_type": "gbdt",
|
|
"num_leaves": 15,
|
|
"learning_rate": 0.03,
|
|
"feature_fraction": 0.7,
|
|
"bagging_fraction": 0.7,
|
|
"bagging_freq": 5,
|
|
"min_data_in_leaf": 50,
|
|
"min_gain_to_split": 0.05,
|
|
"lambda_l1": 0.5,
|
|
"lambda_l2": 1.0,
|
|
"max_depth": 4,
|
|
"verbose": -1,
|
|
"random_state": 42,
|
|
}
|
|
|
|
|
|
class ProbabilityCalibrator:
|
|
"""
|
|
Post-hoc probability calibration using Platt scaling + isotonic regression.
|
|
|
|
Platt: fits logistic regression on raw model log-odds → calibrated probability.
|
|
Works well when raw scores follow a sigmoidal miscalibration pattern.
|
|
Isotonic: non-parametric, fits step-wise monotonic function.
|
|
Better for non-sigmoidal patterns but needs more data.
|
|
|
|
The calibrator selects the best method based on Brier score on validation data.
|
|
"""
|
|
|
|
def __init__(self, min_obs_isotonic: int = 100):
|
|
self.min_obs_isotonic = min_obs_isotonic
|
|
self.platt_model: Optional[LogisticRegression] = None
|
|
self.iso_model: Optional[IsotonicRegression] = None
|
|
self.method: Optional[str] = None # "platt", "isotonic", or "none"
|
|
self.fitted: bool = False
|
|
|
|
def fit(self, raw_scores: np.ndarray, y_true: np.ndarray):
|
|
"""
|
|
Fit calibration on validation data.
|
|
|
|
Parameters
|
|
----------
|
|
raw_scores : np.ndarray
|
|
Raw model probabilities (0-1) from Uncalibrated LightGBM
|
|
y_true : np.ndarray
|
|
Binary ground truth labels
|
|
"""
|
|
if len(raw_scores) < 10:
|
|
self.method = "none"
|
|
self.fitted = True
|
|
return
|
|
|
|
raw_scores = np.clip(raw_scores, 0.001, 0.999).reshape(-1, 1)
|
|
y_true = np.asarray(y_true).ravel()
|
|
|
|
from sklearn.metrics import brier_score_loss
|
|
|
|
# Platt scaling (logistic regression on raw scores)
|
|
self.platt_model = LogisticRegression(C=1.0, solver="lbfgs")
|
|
self.platt_model.fit(raw_scores, y_true)
|
|
platt_proba = self.platt_model.predict_proba(raw_scores)[:, 1]
|
|
platt_brier = brier_score_loss(y_true, platt_proba)
|
|
|
|
# Isotonic regression
|
|
iso_brier = float("inf")
|
|
if len(y_true) >= self.min_obs_isotonic and IsotonicRegression is not None:
|
|
try:
|
|
self.iso_model = IsotonicRegression(
|
|
y_min=0.001, y_max=0.999, out_of_bounds="clip"
|
|
)
|
|
self.iso_model.fit(raw_scores.ravel(), y_true)
|
|
iso_proba = self.iso_model.predict(raw_scores.ravel())
|
|
iso_brier = brier_score_loss(y_true, iso_proba)
|
|
except Exception:
|
|
self.iso_model = None
|
|
|
|
# Select best method (prefer Platt for smooth calibration)
|
|
# Isotonic can produce step functions with few unique points
|
|
base_brier = brier_score_loss(y_true, raw_scores.ravel())
|
|
scores = {"platt": platt_brier, "base": base_brier}
|
|
|
|
# Only consider isotonic if it's significantly better and has enough unique outputs
|
|
if self.iso_model is not None and iso_brier < platt_brier * 0.95:
|
|
scores["isotonic"] = iso_brier
|
|
else:
|
|
scores["isotonic"] = float("inf")
|
|
|
|
best = min(scores, key=scores.get)
|
|
|
|
if best == "isotonic" and self.iso_model is not None:
|
|
self.method = "isotonic"
|
|
elif best == "platt" and self.platt_model is not None:
|
|
self.method = "platt"
|
|
else:
|
|
self.method = "none" # Raw scores are already best
|
|
|
|
self.fitted = True
|
|
print(f" Calibration: {self.method} (platt_brier={platt_brier:.4f}, "
|
|
f"iso_brier={iso_brier:.4f}, raw_brier={base_brier:.4f})")
|
|
|
|
def calibrate(self, raw_scores: np.ndarray) -> np.ndarray:
|
|
"""Apply fitted calibration to raw scores (0-1)."""
|
|
if not self.fitted or self.method == "none":
|
|
raw = np.clip(raw_scores, 0.01, 0.99)
|
|
return np.clip(raw, 0.01, 0.99)
|
|
|
|
raw = np.atleast_1d(raw_scores)
|
|
raw_clipped = np.clip(raw, 0.001, 0.999)
|
|
|
|
if self.method == "platt" and self.platt_model is not None:
|
|
cal = self.platt_model.predict_proba(raw_clipped.reshape(-1, 1))[:, 1]
|
|
elif self.method == "isotonic" and self.iso_model is not None:
|
|
cal = self.iso_model.predict(raw_clipped.ravel())
|
|
else:
|
|
cal = raw_clipped.ravel()
|
|
|
|
# Gentle blending toward 0.5 for extreme probabilities
|
|
# Only blend when raw is very extreme (>0.95 or <0.05)
|
|
extremes = np.abs(raw_clipped.ravel() - 0.5)
|
|
blend = np.clip((extremes - 0.4) / 0.1, 0, 0.3)
|
|
cal_smoothed = cal * (1 - blend) + 0.5 * blend
|
|
|
|
return np.clip(cal_smoothed, 0.01, 0.99)
|
|
|
|
def save(self, path: str):
|
|
"""Save calibration params."""
|
|
data = {
|
|
"method": self.method,
|
|
"platt": pickle.dumps(self.platt_model) if self.platt_model else None,
|
|
"iso": pickle.dumps(self.iso_model) if self.iso_model else None,
|
|
}
|
|
with open(path, "wb") as f:
|
|
pickle.dump(data, f)
|
|
|
|
def load(self, path: str):
|
|
"""Load calibration params."""
|
|
with open(path, "rb") as f:
|
|
data = pickle.load(f)
|
|
self.method = data.get("method", "none")
|
|
if data.get("platt"):
|
|
self.platt_model = pickle.loads(data["platt"])
|
|
if data.get("iso"):
|
|
self.iso_model = pickle.loads(data["iso"])
|
|
self.fitted = True
|
|
|
|
|
|
class WeatherModel:
|
|
"""Probability model for a single weather target.
|
|
|
|
Two model modes:
|
|
- 'lgb': LightGBM gradient boosting (for real ERA5 data)
|
|
- 'lr': Logistic regression (for synthetic/bootstrap data, prevents overfitting)
|
|
"""
|
|
|
|
def __init__(self, target_name: str, mode: str = "lgb"):
|
|
if target_name not in TARGET_DEFINITIONS:
|
|
raise ValueError(f"Unknown target: {target_name}")
|
|
self.target_name = target_name
|
|
self.target_def = TARGET_DEFINITIONS[target_name]
|
|
self.model: Optional[lgb.Booster] = None
|
|
self.lr_model = None # LogisticRegression for 'lr' mode
|
|
self.mode = mode
|
|
self.feature_importance: Dict[str, float] = {}
|
|
self.calibrator = ProbabilityCalibrator()
|
|
self._trained = False
|
|
|
|
def train(
|
|
self,
|
|
X_train: np.ndarray,
|
|
y_train: np.ndarray,
|
|
X_val: Optional[np.ndarray] = None,
|
|
y_val: Optional[np.ndarray] = None,
|
|
params: Optional[Dict] = None,
|
|
early_stopping_rounds: int = 50,
|
|
verbose: bool = True,
|
|
):
|
|
"""Train model + calibrate."""
|
|
if self.mode == "lr":
|
|
self._train_lr(X_train, y_train, X_val, y_val)
|
|
else:
|
|
self._train_lgb(X_train, y_train, X_val, y_val, params, early_stopping_rounds, verbose)
|
|
|
|
self._trained = True
|
|
if X_val is not None and y_val is not None:
|
|
self.calibrator.fit(self.predict_raw(X_val), y_val)
|
|
elif X_train is not None and y_train is not None:
|
|
self.calibrator.fit(self.predict_raw(X_train), y_train)
|
|
|
|
def _train_lr(self, X_train, y_train, X_val, y_val):
|
|
"""Train logistic regression model."""
|
|
if LogisticRegression is None:
|
|
raise ImportError("scikit-learn not installed")
|
|
from sklearn.preprocessing import StandardScaler
|
|
self.scaler = StandardScaler()
|
|
X_train_scaled = self.scaler.fit_transform(X_train)
|
|
self.lr_model = LogisticRegression(
|
|
C=0.1, # Strong L2 regularization
|
|
solver="lbfgs",
|
|
max_iter=2000,
|
|
class_weight="balanced",
|
|
)
|
|
self.lr_model.fit(X_train_scaled, y_train)
|
|
self._trained = True
|
|
|
|
def _train_lgb(self, X_train, y_train, X_val, y_val, params, early_stopping_rounds, verbose):
|
|
"""Train LightGBM model."""
|
|
if lgb is None:
|
|
raise ImportError("lightgbm not installed")
|
|
train_params = {**LGBM_PARAMS, **(params or {})}
|
|
dtrain = lgb.Dataset(X_train, label=y_train)
|
|
if X_val is not None and y_val is not None:
|
|
dval = lgb.Dataset(X_val, label=y_val, reference=dtrain)
|
|
valid_sets, valid_names = [dtrain, dval], ["train", "valid"]
|
|
else:
|
|
valid_sets, valid_names = None, None
|
|
self.model = lgb.train(
|
|
train_params, dtrain, num_boost_round=500,
|
|
valid_sets=valid_sets, valid_names=valid_names,
|
|
callbacks=[lgb.early_stopping(early_stopping_rounds),
|
|
lgb.log_evaluation(period=50 if verbose else 0)]
|
|
if X_val is not None else None,
|
|
)
|
|
self._compute_feature_importance()
|
|
|
|
def predict_raw(self, X: np.ndarray) -> np.ndarray:
|
|
"""Raw probability (0-1) before calibration."""
|
|
if not self._trained:
|
|
raise RuntimeError("Model not trained")
|
|
if self.mode == "lr" and self.lr_model is not None:
|
|
X_scaled = self.scaler.transform(X)
|
|
return self.lr_model.predict_proba(X_scaled)[:, 1]
|
|
elif self.model is not None:
|
|
return self.model.predict(X)
|
|
else:
|
|
return np.full(len(X), 0.5)
|
|
|
|
def predict_proba(self, X: np.ndarray) -> np.ndarray:
|
|
"""Calibrated probability (0-100)."""
|
|
raw = self.predict_raw(X)
|
|
cal = self.calibrator.calibrate(raw)
|
|
return cal * 100.0
|
|
|
|
def predict(self, X: np.ndarray, threshold: float = 50.0) -> np.ndarray:
|
|
return (self.predict_proba(X) >= threshold).astype(int)
|
|
|
|
def evaluate(self, X: np.ndarray, y: np.ndarray) -> Dict[str, float]:
|
|
proba = self.predict_proba(X) / 100.0
|
|
pred = (proba >= 0.5).astype(int)
|
|
from sklearn.metrics import accuracy_score, brier_score_loss, roc_auc_score, log_loss
|
|
return {
|
|
"accuracy": float(accuracy_score(y, pred)),
|
|
"brier_score": float(brier_score_loss(y, proba)),
|
|
"roc_auc": float(roc_auc_score(y, proba)) if len(np.unique(y)) > 1 else 0.5,
|
|
"log_loss": float(log_loss(y, proba)),
|
|
"n_samples": len(y),
|
|
"p_yes_actual": float(y.mean() * 100),
|
|
"p_yes_predicted": float(proba.mean() * 100),
|
|
"calibration_method": self.calibrator.method,
|
|
}
|
|
|
|
def _compute_feature_importance(self):
|
|
if self.model is None:
|
|
return
|
|
gain = self.model.feature_importance(importance_type="gain")
|
|
names = self.model.feature_name()
|
|
self.feature_importance = dict(sorted(zip(names, gain), key=lambda x: x[1], reverse=True))
|
|
|
|
def top_features(self, n: int = 15) -> Dict[str, float]:
|
|
items = sorted(self.feature_importance.items(), key=lambda x: x[1], reverse=True)
|
|
return dict(items[:n])
|
|
|
|
def save(self, path: Optional[str] = None):
|
|
MODEL_DIR.mkdir(parents=True, exist_ok=True)
|
|
p = path or (MODEL_DIR / f"{self.target_name}.lgb")
|
|
if self.model:
|
|
self.model.save_model(str(p))
|
|
elif not os.path.exists(p):
|
|
# Create marker for LR models
|
|
with open(p, "w") as f:
|
|
f.write("lr")
|
|
|
|
meta = {
|
|
"target_name": self.target_name,
|
|
"target_definition": self.target_def,
|
|
"mode": self.mode,
|
|
"feature_importance": self.feature_importance,
|
|
"trained": self._trained,
|
|
"calibration_method": self.calibrator.method,
|
|
}
|
|
with open(str(p).replace(".lgb", "_meta.json"), "w") as f:
|
|
json.dump(meta, f, indent=2, default=str)
|
|
self.calibrator.save(str(p).replace(".lgb", "_cal.pkl"))
|
|
if self.lr_model is not None:
|
|
import pickle
|
|
with open(str(p).replace(".lgb", "_lr.pkl"), "wb") as f:
|
|
pickle.dump({"model": self.lr_model, "scaler": self.scaler}, f)
|
|
|
|
def load(self, path: Optional[str] = None):
|
|
p = path or (MODEL_DIR / f"{self.target_name}.lgb")
|
|
if not os.path.exists(p):
|
|
raise FileNotFoundError(f"Model not found: {p}")
|
|
meta_path = str(p).replace(".lgb", "_meta.json")
|
|
if os.path.exists(meta_path):
|
|
with open(meta_path) as f:
|
|
meta = json.load(f)
|
|
self.mode = meta.get("mode", "lgb")
|
|
self.feature_importance = meta.get("feature_importance", {})
|
|
|
|
if self.mode == "lr":
|
|
import pickle
|
|
lr_path = str(p).replace(".lgb", "_lr.pkl")
|
|
if os.path.exists(lr_path):
|
|
with open(lr_path, "rb") as f:
|
|
data = pickle.load(f)
|
|
self.lr_model = data["model"]
|
|
self.scaler = data["scaler"]
|
|
else:
|
|
if lgb is None:
|
|
raise ImportError("lightgbm not installed")
|
|
self.model = lgb.Booster(model_file=str(p))
|
|
self._trained = True
|
|
cal_path = str(p).replace(".lgb", "_cal.pkl")
|
|
if os.path.exists(cal_path):
|
|
self.calibrator.load(cal_path)
|
|
|
|
@staticmethod
|
|
def build_target(df: pd.DataFrame, variable: str, threshold: float, op: str = "gt") -> np.ndarray:
|
|
values = df[variable].values
|
|
if op == "gt":
|
|
return (values > threshold).astype(int)
|
|
elif op == "ge":
|
|
return (values >= threshold).astype(int)
|
|
elif op == "lt":
|
|
return (values < threshold).astype(int)
|
|
elif op == "le":
|
|
return (values <= threshold).astype(int)
|
|
else:
|
|
raise ValueError(f"Unknown operator: {op}")
|
|
|
|
|
|
class ModelEnsemble:
|
|
"""Manage multiple WeatherModel instances."""
|
|
|
|
def __init__(self):
|
|
self.models: Dict[str, WeatherModel] = {}
|
|
|
|
def load_all(self):
|
|
MODEL_DIR.mkdir(parents=True, exist_ok=True)
|
|
for target in TARGET_DEFINITIONS:
|
|
model_path = MODEL_DIR / f"{target}.lgb"
|
|
if model_path.exists():
|
|
model = WeatherModel(target)
|
|
model.load(str(model_path))
|
|
self.models[target] = model
|
|
return self.models
|
|
|
|
def predict_all(self, X: np.ndarray) -> Dict[str, float]:
|
|
if X.ndim == 1:
|
|
X = X.reshape(1, -1)
|
|
return {name: float(model.predict_proba(X)[0]) for name, model in self.models.items()}
|
|
|
|
def predict(self, target: str, X: np.ndarray) -> float:
|
|
if target not in self.models:
|
|
raise KeyError(f"Model '{target}' not loaded.")
|
|
return float(self.models[target].predict_proba(X)[0])
|
|
|
|
def has(self, target: str) -> bool:
|
|
return target in self.models
|
|
|
|
@property
|
|
def available_targets(self) -> List[str]:
|
|
return list(self.models.keys())
|