Machine Learning: 本番環境のモニタリングとモデルドリフト — デプロイはゴールではなくスタートライン

def check_drift(reference_data, current_data, threshold=0.2):
    psi = calculate_psi(reference_data, current_data)
    if psi > threshold:
        alert("Data drift detected! PSI={:.3f}".format(psi))
        trigger_retraining()

import numpy as np
from datetime import datetime, timedelta

class ModelMonitor:
    def __init__(self):
        self.metrics_history = []

    def log_prediction(self, prediction, features, latency_ms):
        self.metrics_history.append({
            "timestamp": datetime.now(),
            "prediction": prediction,
            "features_mean": np.mean(features),
            "latency_ms": latency_ms,
        })

    def check_health(self, last_n=1000):
        if len(self.metrics_history) < last_n:
            return "Insufficient data"

        recent = self.metrics_history[-last_n:]
        predictions = [m["prediction"] for m in recent]
        latencies = [m["latency_ms"] for m in recent]

        report = {
            "prediction_mean": np.mean(predictions),
            "prediction_std": np.std(predictions),
            "p95_latency_ms": np.percentile(latencies, 95),
            "sample_count": len(recent),
        }
        return report

monitor = ModelMonitor()
for _ in range(100):
    monitor.log_prediction(prediction=np.random.normal(200, 50),
                           features=np.random.randn(5),
                           latency_ms=np.random.uniform(10, 30))
print(monitor.check_health())

# Functions defined successfully

import numpy as np

def calculate_psi(reference, current, n_bins=10):
    """Calculate Population Stability Index.

    PSI < 0.1: No significant change
    PSI 0.1-0.2: Moderate change, investigate
    PSI > 0.2: Significant change, action required
    """
    breakpoints = np.percentile(reference, np.linspace(0, 100, n_bins + 1))
    breakpoints[0] = -np.inf
    breakpoints[-1] = np.inf

    ref_counts = np.histogram(reference, bins=breakpoints)[0]
    cur_counts = np.histogram(current, bins=breakpoints)[0]

    ref_pct = ref_counts / len(reference)
    cur_pct = cur_counts / len(current)

    # Avoid log(0)
    ref_pct = np.clip(ref_pct, 1e-6, None)
    cur_pct = np.clip(cur_pct, 1e-6, None)

    psi = np.sum((cur_pct - ref_pct) * np.log(cur_pct / ref_pct))
    return psi

# Simulate: normal vs drifted data
rng = np.random.default_rng(42)
reference = rng.normal(100, 20, 10000)

# No drift
current_no_drift = rng.normal(100, 20, 5000)
psi_no_drift = calculate_psi(reference, current_no_drift)

# Moderate drift (mean shifted by 10)
current_moderate = rng.normal(110, 20, 5000)
psi_moderate = calculate_psi(reference, current_moderate)

# Severe drift (mean shifted by 25, Double-11 effect)
current_severe = rng.normal(125, 30, 5000)
psi_severe = calculate_psi(reference, current_severe)

print(f"No drift:    PSI = {psi_no_drift:.3f} (OK)")
print(f"Moderate:    PSI = {psi_moderate:.3f} (Investigate)")
print(f"Severe:      PSI = {psi_severe:.3f} (Action Required!)")

# Functions defined successfully

from scipy import stats
import numpy as np

def calculate_js_divergence(reference, current, n_bins=50):
    """Jensen-Shannon Divergence (symmetric, bounded [0, 1])."""
    bins = np.linspace(min(reference.min(), current.min()),
                       max(reference.max(), current.max()), n_bins + 1)
    p = np.histogram(reference, bins=bins, density=True)[0]
    q = np.histogram(current, bins=bins, density=True)[0]
    m = (p + q) / 2
    js = 0.5 * stats.entropy(p, m) + 0.5 * stats.entropy(q, m)
    return js

rng = np.random.default_rng(42)
ref = rng.normal(100, 20, 10000)
cur_drifted = rng.normal(115, 25, 5000)

# KS test
ks_stat, ks_pvalue = stats.ks_2samp(ref, cur_drifted)

# PSI
psi = calculate_psi(ref, cur_drifted)

# JS divergence
js = calculate_js_divergence(ref, cur_drifted)

print(f"KS statistic: {ks_stat:.4f}, p-value: {ks_pvalue:.6f}")
print(f"PSI:          {psi:.4f}")
print(f"JS Divergence: {js:.4f}")

# Functions defined successfully

import numpy as np

rng = np.random.default_rng(42)

# Normal period: ad_spend → revenue relationship
n_normal = 10000
ad_spend_normal = rng.uniform(10, 100, n_normal)
revenue_normal = 50 + 0.8 * ad_spend_normal + rng.normal(0, 10, n_normal)

# Double-11 period: behavior completely changed
n_promo = 5000
ad_spend_promo = rng.uniform(50, 200, n_promo)  # Higher ad spend
revenue_promo = 200 + 1.5 * ad_spend_promo + rng.normal(0, 30, n_promo)  # Different coefficients!

# Detect drift in ad_spend distribution
psi = calculate_psi(ad_spend_normal, ad_spend_promo)
ks_stat, ks_p = stats.ks_2samp(ad_spend_normal, ad_spend_promo)

print("=== Data Drift Detection ===")
print(f"Ad spend PSI: {psi:.3f} {'⚠️ DRIFT' if psi > 0.2 else 'OK'}")
print(f"Ad spend KS:  {ks_stat:.3f} (p={ks_p:.6f}) {'⚠️ DRIFT' if ks_p < 0.05 else 'OK'}")

# Detect concept drift: same X, different Y|X relationship
# Use a model trained on normal data to predict promo data
from sklearn.linear_model import LinearRegression
model_normal = LinearRegression()
model_normal.fit(ad_spend_normal.reshape(-1, 1), revenue_normal)

pred_promo = model_normal.predict(ad_spend_promo.reshape(-1, 1))
mape_normal_on_normal = np.mean(np.abs((revenue_normal - model_normal.predict(ad_spend_normal.reshape(-1, 1))) / revenue_normal)) * 100
mape_normal_on_promo = np.mean(np.abs((revenue_promo - pred_promo) / revenue_promo)) * 100

print(f"\n=== Concept Drift Detection ===")
print(f"MAPE (normal→normal): {mape_normal_on_normal:.1f}%")
print(f"MAPE (normal→promo):  {mape_normal_on_promo:.1f}%")
print(f"Performance degradation: {mape_normal_on_promo - mape_normal_on_normal:.1f}% ⚠️ CONCEPT DRIFT")

=== Data Drift Detection ===

import numpy as np
from sklearn.linear_model import LinearRegression
from sklearn.metrics import mean_absolute_percentage_error

class AutoRetrainPipeline:
    def __init__(self, drift_threshold_psi=0.2, mape_threshold=0.15):
        self.drift_threshold_psi = drift_threshold_psi
        self.mape_threshold = mape_threshold
        self.model = None
        self.reference_data = None

    def train(self, X, y):
        self.model = LinearRegression()
        self.model.fit(X, y)
        self.reference_data = X.copy()

    def check_drift(self, X_current):
        drifts = []
        for i in range(X_current.shape[1]):
            psi = calculate_psi(self.reference_data[:, i], X_current[:, i])
            drifts.append({"feature": f"feat_{i}", "psi": psi, "drift": psi > self.drift_threshold_psi})
        return drifts

    def check_performance(self, X, y):
        y_pred = self.model.predict(X)
        mape = mean_absolute_percentage_error(y, y_pred)
        return mape < self.mape_threshold, mape

    def auto_retrain_if_needed(self, X_current, y_current):
        # Step 1: Check drift
        drifts = self.check_drift(X_current)
        any_drift = any(d["drift"] for d in drifts)

        # Step 2: Check performance
        performance_ok, mape = self.check_performance(X_current, y_current)

        # Step 3: Decision
        if not any_drift and performance_ok:
            return "OK - No retraining needed"

        if any_drift and not performance_ok:
            # Retrain with recent data
            self.model.fit(X_current, y_current)
            self.reference_data = X_current.copy()
            new_mape = mean_absolute_percentage_error(y_current, self.model.predict(X_current))
            return f"RETRAINED - PSI drift detected, MAPE {mape:.1%} → {new_mape:.1%}"

        if any_drift:
            return f"DRIFT WARNING - PSI drift but performance still OK (MAPE {mape:.1%})"

        return f"PERFORMANCE WARNING - No drift but MAPE degraded to {mape:.1%}"

# Test
pipeline = AutoRetrainPipeline()
rng = np.random.default_rng(42)
X_train = rng.uniform(10, 100, (5000, 3))
y_train = 50 + 0.8 * X_train[:, 0] + rng.normal(0, 5, 5000)
pipeline.train(X_train, y_train)

# Simulate Double-11 data
X_promo = rng.uniform(50, 200, (2000, 3))
y_promo = 200 + 1.5 * X_promo[:, 0] + rng.normal(0, 10, 2000)

result = pipeline.auto_retrain_if_needed(X_promo, y_promo)
print(result)

# Functions defined successfully