05 · Classical ML in Production

The gap between a working notebook and a deployed, maintainable model is where most ML projects die. This file covers the mechanics of taking a trained sklearn model and making it: reproducible, serializable, serveable, and monitored. None of this is glamorous work. All of it is necessary work. A model nobody can use is a research artifact, not an engineering product.

⚠️ What Most People Get Wrong: Treating serialization as the finish line. Pickling a model and handing it to an API is day one, not the end. The real production questions are: Does it handle malformed inputs? Will it degrade gracefully when the input distribution shifts? Can you retrain it without re-architecting the service? This file answers those questions.

---

1. The sklearn Pipeline — Your Serialization Contract

Before you can serialize anything useful, your preprocessing and model must be wrapped in a single sklearn Pipeline. This is not optional. A pipeline serializes the entire transformation chain — imputers, scalers, encoders, the model — as one atomic unit. Serializing only the model and not the preprocessors is the most common cause of train-serve skew.

import pandas as pd
import numpy as np
from sklearn.pipeline import Pipeline
from sklearn.compose import ColumnTransformer
from sklearn.preprocessing import StandardScaler, OneHotEncoder
from sklearn.impute import SimpleImputer
from sklearn.ensemble import GradientBoostingClassifier
from sklearn.model_selection import train_test_split
from sklearn.metrics import roc_auc_score

# --- Build a production-ready pipeline ---

# Column definitions (define these once, use everywhere)
NUMERIC_FEATURES = ['age', 'income', 'tenure_days', 'num_transactions', 'avg_order_value']
CATEGORICAL_FEATURES = ['country', 'product_category', 'acquisition_channel']
TARGET = 'churned'

# Preprocessing
numeric_transformer = Pipeline([
    ('imputer', SimpleImputer(strategy='median')),
    ('scaler', StandardScaler())
])

categorical_transformer = Pipeline([
    ('imputer', SimpleImputer(strategy='constant', fill_value='MISSING')),
    ('encoder', OneHotEncoder(handle_unknown='ignore', sparse_output=False))
])

preprocessor = ColumnTransformer([
    ('num', numeric_transformer, NUMERIC_FEATURES),
    ('cat', categorical_transformer, CATEGORICAL_FEATURES)
])

# Full pipeline
pipeline = Pipeline([
    ('preprocessor', preprocessor),
    ('classifier', GradientBoostingClassifier(
        n_estimators=200,
        learning_rate=0.05,
        max_depth=4,
        subsample=0.8,
        random_state=42
    ))
])

# Train
pipeline.fit(X_train, y_train)

# Evaluate
y_prob = pipeline.predict_proba(X_test)[:, 1]
print(f"Test AUC: {roc_auc_score(y_test, y_prob):.4f}")

# The pipeline handles all input types correctly at inference time:
# single row prediction (returns probability)
sample = X_test.iloc[[0]]  # Keep as DataFrame, not Series
prob = pipeline.predict_proba(sample)[0, 1]
print(f"Churn probability for sample: {prob:.3f}")

---

2. Model Serialization: joblib vs pickle

joblib (Preferred for sklearn)

import joblib
import os
from datetime import datetime

MODEL_DIR = "models/"
os.makedirs(MODEL_DIR, exist_ok=True)

# Save with version tag
version = datetime.now().strftime("%Y%m%d_%H%M%S")
model_path = f"{MODEL_DIR}/churn_model_{version}.joblib"

joblib.dump(pipeline, model_path, compress=3)  # compress=3 is a good balance
file_size_mb = os.path.getsize(model_path) / (1024 * 1024)
print(f"Model saved to {model_path} ({file_size_mb:.2f} MB)")

# Load
loaded_pipeline = joblib.load(model_path)
# Verify round-trip
assert roc_auc_score(y_test, loaded_pipeline.predict_proba(X_test)[:, 1]) == \
       roc_auc_score(y_test, pipeline.predict_proba(X_test)[:, 1])
print("Serialization round-trip verified ✅")

Why Not pickle Directly?

import pickle

# pickle works but has two problems for sklearn:
# 1. Large models (forests): joblib uses memory-mapped numpy arrays — much faster
# 2. Python version sensitivity: pickled sklearn objects can break across Python minor versions

# Benchmark (approximate for a RF with 200 trees, n=50k, p=100):
# pickle: ~2.1s save, ~1.8s load, ~45MB
# joblib (compress=0): ~0.3s save, ~0.2s load, ~45MB (fast I/O)
# joblib (compress=3): ~0.4s save, ~0.3s load, ~12MB (good balance)
# joblib (compress=9): ~2.0s save, ~0.3s load, ~9MB (diminishing returns)

# Rule: use joblib for sklearn. Use pickle only for non-sklearn objects.

Model Metadata Contract

Always save metadata alongside the model. A model file with no metadata is undeployable in six months.

import json

metadata = {
    "model_version": version,
    "model_type": "GradientBoostingClassifier",
    "training_date": datetime.now().isoformat(),
    "python_version": "3.10.12",
    "sklearn_version": "1.3.0",
    "training_data": {
        "n_samples": len(X_train),
        "n_features": len(NUMERIC_FEATURES) + len(CATEGORICAL_FEATURES),
        "target": TARGET,
        "class_balance": float(y_train.mean())
    },
    "evaluation": {
        "test_auc": float(roc_auc_score(y_test, y_prob)),
        "test_n": len(y_test)
    },
    "features": {
        "numeric": NUMERIC_FEATURES,
        "categorical": CATEGORICAL_FEATURES
    },
    "hyperparameters": pipeline.named_steps['classifier'].get_params()
}

meta_path = model_path.replace('.joblib', '_metadata.json')
with open(meta_path, 'w') as f:
    json.dump(metadata, f, indent=2)

print(f"Metadata saved to {meta_path}")

---

3. Serving a Model as a REST API (FastAPI)

FastAPI is the correct choice for serving sklearn models in 2024+. Flask is viable but has no async support or automatic validation. The performance difference for sklearn (CPU-bound, synchronous) is small, but the developer experience of FastAPI is significantly better.

# File: serve_model.py
# Run with: uvicorn serve_model:app --host 0.0.0.0 --port 8000

from fastapi import FastAPI, HTTPException
from pydantic import BaseModel, Field, validator
from typing import Optional
import joblib
import pandas as pd
import numpy as np
import logging
from datetime import datetime
import time

# Configure logging
logging.basicConfig(level=logging.INFO)
logger = logging.getLogger(__name__)

app = FastAPI(
    title="Churn Prediction API",
    description="Predicts customer churn probability using a GBM model",
    version="1.0.0"
)

# Load model at startup (not per-request — that would be catastrophic)
MODEL_PATH = "models/churn_model_latest.joblib"
pipeline = joblib.load(MODEL_PATH)
logger.info(f"Model loaded from {MODEL_PATH}")


# --- Request/Response Schemas ---

class ChurnPredictionRequest(BaseModel):
    age: float = Field(..., ge=0, le=120, description="Customer age in years")
    income: float = Field(..., ge=0, description="Annual income in USD")
    tenure_days: int = Field(..., ge=0, description="Days since customer joined")
    num_transactions: int = Field(..., ge=0)
    avg_order_value: float = Field(..., ge=0)
    country: str = Field(..., min_length=2, max_length=50)
    product_category: str
    acquisition_channel: str
    
    class Config:
        json_schema_extra = {
            "example": {
                "age": 35.0,
                "income": 75000.0,
                "tenure_days": 365,
                "num_transactions": 12,
                "avg_order_value": 87.50,
                "country": "US",
                "product_category": "electronics",
                "acquisition_channel": "organic_search"
            }
        }


class ChurnPredictionResponse(BaseModel):
    churn_probability: float
    churn_predicted: bool
    confidence: str  # 'HIGH' | 'MEDIUM' | 'LOW'
    model_version: str
    inference_time_ms: float


@app.get("/health")
def health_check():
    """Liveness probe — always returns 200 if the service is running."""
    return {"status": "ok", "timestamp": datetime.now().isoformat()}


@app.get("/model/info")
def model_info():
    """Returns metadata about the loaded model."""
    try:
        import json
        meta_path = MODEL_PATH.replace('.joblib', '_metadata.json')
        with open(meta_path) as f:
            return json.load(f)
    except FileNotFoundError:
        return {"error": "Metadata file not found"}


@app.post("/predict", response_model=ChurnPredictionResponse)
def predict_churn(request: ChurnPredictionRequest):
    """
    Returns churn probability for a single customer.
    
    Confidence levels:
    - HIGH: probability < 0.2 or > 0.8
    - MEDIUM: probability 0.2-0.35 or 0.65-0.8
    - LOW: probability 0.35-0.65 (model is uncertain)
    """
    start_time = time.time()
    
    try:
        # Convert to DataFrame (preserves feature names for pipeline)
        input_df = pd.DataFrame([request.dict()])
        
        # Predict
        prob = float(pipeline.predict_proba(input_df)[0, 1])
        predicted = prob >= 0.5
        
        # Confidence tier
        if prob < 0.2 or prob > 0.8:
            confidence = "HIGH"
        elif prob < 0.35 or prob > 0.65:
            confidence = "MEDIUM"
        else:
            confidence = "LOW"
        
        inference_ms = (time.time() - start_time) * 1000
        logger.info(f"Prediction: prob={prob:.4f}, latency={inference_ms:.2f}ms")
        
        return ChurnPredictionResponse(
            churn_probability=round(prob, 4),
            churn_predicted=predicted,
            confidence=confidence,
            model_version="1.0.0",
            inference_time_ms=round(inference_ms, 2)
        )
    
    except Exception as e:
        logger.error(f"Prediction error: {e}")
        raise HTTPException(status_code=500, detail=f"Prediction failed: {str(e)}")


@app.post("/predict/batch")
def predict_batch(requests: list[ChurnPredictionRequest]):
    """Batch prediction — more efficient than N individual calls."""
    if len(requests) > 1000:
        raise HTTPException(status_code=400, detail="Batch size limit: 1000")
    
    start_time = time.time()
    input_df = pd.DataFrame([r.dict() for r in requests])
    probs = pipeline.predict_proba(input_df)[:, 1].tolist()
    inference_ms = (time.time() - start_time) * 1000
    
    return {
        "predictions": [{"churn_probability": round(p, 4), "churn_predicted": p >= 0.5}
                        for p in probs],
        "batch_size": len(requests),
        "total_inference_time_ms": round(inference_ms, 2),
        "avg_per_sample_ms": round(inference_ms / len(requests), 2)
    }

Running the API

# Install dependencies
pip install fastapi uvicorn[standard] joblib pandas scikit-learn

# Run development server (auto-reload)
uvicorn serve_model:app --reload --host 0.0.0.0 --port 8000

# Test it
curl -X POST "http://localhost:8000/predict" \
  -H "Content-Type: application/json" \
  -d '{"age": 35, "income": 75000, "tenure_days": 365, "num_transactions": 12,
       "avg_order_value": 87.5, "country": "US", "product_category": "electronics",
       "acquisition_channel": "organic_search"}'

# Expected response (~5ms latency for GBM, single prediction):
# {"churn_probability": 0.1823, "churn_predicted": false, "confidence": "HIGH", ...}

# Interactive docs (auto-generated by FastAPI):
# http://localhost:8000/docs

Performance Characteristics (Reference Numbers)

For a GradientBoostingClassifier (n_estimators=200, n_features=50, tabular data):

• Single prediction latency: ~2-8ms (CPU, no GPU needed for sklearn)

• Batch of 100: ~15-25ms total (~0.15-0.25ms per sample — batching helps)

• Throughput: ~500-1000 RPS on a single CPU core for simple GBM

• Random Forest (200 trees): ~5-15ms single prediction (more trees = slower)

• XGBoost: ~1-3ms single prediction (typically fastest sklearn-equivalent)

---

4. Model Versioning Basics

models/
├── churn_model_20241015_143022.joblib
├── churn_model_20241015_143022_metadata.json
├── churn_model_20241201_091155.joblib
├── churn_model_20241201_091155_metadata.json
└── churn_model_latest.joblib  ← symlink to current production model

import os
import shutil

def promote_to_production(model_path: str, models_dir: str = "models/"):
    """
    Promote a specific model version to production.
    Creates/updates the 'latest' symlink.
    """
    latest_path = os.path.join(models_dir, "churn_model_latest.joblib")
    
    # Backup current production model
    if os.path.exists(latest_path):
        backup_path = latest_path.replace('latest', 'prev_production')
        shutil.copy2(latest_path, backup_path)
        print(f"Previous production backed up to {backup_path}")
    
    # Copy new version to latest
    shutil.copy2(model_path, latest_path)
    print(f"Promoted {model_path} to production ({latest_path})")
    
    # Verify
    loaded = joblib.load(latest_path)
    print(f"Verification: {type(loaded).__name__} loaded successfully")

# Usage
promote_to_production("models/churn_model_20241201_091155.joblib")

---

5. Monitoring for Prediction Drift

A model that was 85% accurate at launch can degrade silently as the real-world data distribution shifts. Without monitoring, you find out when someone complains, not when the shift starts.

The Two Types of Drift

Type	What It Is	Example	How to Detect
Data Drift (Covariate Shift)	P(X) changes, P(Y|X) stays same	Average customer age shifts after new campaign	Statistical tests on feature distributions
Concept Drift	P(Y|X) changes	Fraud patterns change after attackers adapt	Monitor prediction accuracy on labeled data

Implementing a Basic Drift Monitor

import numpy as np
import pandas as pd
from scipy import stats

class PredictionDriftMonitor:
    """
    Monitors feature distributions and prediction score distributions
    for statistically significant drift vs. a reference window.
    """
    
    def __init__(self, reference_df: pd.DataFrame, 
                 reference_probs: np.ndarray,
                 alpha: float = 0.05):
        """
        reference_df: training/recent-production feature data
        reference_probs: model output probabilities on reference data
        alpha: significance level for drift detection
        """
        self.reference_df = reference_df
        self.reference_probs = reference_probs
        self.alpha = alpha
        self.numeric_cols = reference_df.select_dtypes(include=np.number).columns.tolist()
    
    def check_feature_drift(self, current_df: pd.DataFrame) -> dict:
        """
        Kolmogorov-Smirnov test for each numeric feature.
        KS test is distribution-free: works for any shape.
        """
        results = {}
        drifted_features = []
        
        for col in self.numeric_cols:
            if col not in current_df.columns:
                continue
            
            ref_vals = self.reference_df[col].dropna().values
            cur_vals = current_df[col].dropna().values
            
            ks_stat, p_value = stats.ks_2samp(ref_vals, cur_vals)
            drifted = p_value < self.alpha
            
            results[col] = {
                'ks_statistic': round(ks_stat, 4),
                'p_value': round(p_value, 4),
                'drifted': drifted
            }
            
            if drifted:
                drifted_features.append(col)
        
        return {
            'feature_tests': results,
            'drifted_features': drifted_features,
            'drift_detected': len(drifted_features) > 0,
            'drift_fraction': len(drifted_features) / len(self.numeric_cols)
        }
    
    def check_prediction_drift(self, current_probs: np.ndarray) -> dict:
        """
        Monitor the distribution of model output probabilities.
        Shift here is a strong signal of concept drift.
        """
        ks_stat, p_value = stats.ks_2samp(self.reference_probs, current_probs)
        
        ref_mean = self.reference_probs.mean()
        cur_mean = current_probs.mean()
        mean_shift = cur_mean - ref_mean
        
        psi = self._population_stability_index(self.reference_probs, current_probs)
        
        return {
            'ks_statistic': round(ks_stat, 4),
            'p_value': round(p_value, 4),
            'prediction_drift': p_value < self.alpha,
            'reference_mean_prob': round(ref_mean, 4),
            'current_mean_prob': round(cur_mean, 4),
            'mean_shift': round(mean_shift, 4),
            'psi': round(psi, 4),
            'psi_interpretation': 'stable' if psi < 0.1 else 'moderate change' if psi < 0.25 else 'significant shift'
        }
    
    @staticmethod
    def _population_stability_index(reference: np.ndarray, 
                                     current: np.ndarray, 
                                     bins: int = 10) -> float:
        """
        PSI: industry-standard metric for prediction distribution shift.
        PSI < 0.1: no significant change
        PSI 0.1-0.25: moderate change — investigate
        PSI > 0.25: significant shift — likely retrain needed
        """
        ref_hist, bin_edges = np.histogram(reference, bins=bins, range=(0, 1))
        cur_hist, _ = np.histogram(current, bins=bin_edges)
        
        ref_pct = (ref_hist + 1e-6) / len(reference)
        cur_pct = (cur_hist + 1e-6) / len(current)
        
        psi = np.sum((cur_pct - ref_pct) * np.log(cur_pct / ref_pct))
        return psi


# Usage example
monitor = PredictionDriftMonitor(X_train_df, pipeline.predict_proba(X_train_df)[:, 1])

# Check weekly on new production data
feature_report = monitor.check_feature_drift(X_new_production_df)
prediction_report = monitor.check_prediction_drift(
    pipeline.predict_proba(X_new_production_df)[:, 1]
)

print(f"Drifted features: {feature_report['drifted_features']}")
print(f"PSI: {prediction_report['psi']} ({prediction_report['psi_interpretation']})")

if prediction_report['psi'] > 0.25:
    print("⚠️  PSI > 0.25: Schedule model retraining.")

---

6. The Full Production Checklist

Before handing off a model to production, run through this:

PRE-DEPLOYMENT
□ Pipeline wraps ALL preprocessing (no naked scalers/encoders)
□ Pipeline verified round-trip: fit → save → load → predict matches
□ Model metadata JSON saved alongside .joblib
□ Input validation: schema enforced (Pydantic / pandera)
□ Edge cases tested: missing values, unseen categories, extreme values
□ Latency benchmarked under expected load
□ Memory usage measured (large forests can be 500MB+ in RAM)

API
□ Health check endpoint (/health) returns 200
□ Input validation returns 422 (not 500) on bad input
□ Batch endpoint available for high-throughput use cases
□ Logging: request, latency, prediction score (not features — PII risk)
□ Error handling: 500 errors are logged with context

MONITORING
□ Feature drift monitor running on weekly production windows
□ PSI alert set for > 0.1 (investigate) and > 0.25 (retrain)
□ Ground truth labels collected for periodic accuracy audit
□ Previous model version kept as rollback target

---

7. Bridge to MLOps

This file covers the essentials of production sklearn. The full MLOps discipline — experiment tracking, automated retraining, CI/CD for models, feature stores, A/B testing infrastructure — is Phase 5 of this roadmap. The patterns you’ve built here (pipeline-first, metadata contracts, drift monitoring) are directly portable to MLflow, DVC, and Seldon.

The mental model that carries forward: a model is not done when it achieves good test AUC. A model is done when it can be retrained, versioned, deployed, and monitored without manual intervention.

---

Resources

Resource	Type	Why	Time
FastAPI Docs	Docs	Canonical reference, best async Python API framework	Reference
“Designing ML Systems” — Chip Huyen	Book	Best end-to-end production ML book, covers drift and monitoring	8 hrs
sklearn Pipeline docs	Docs	ColumnTransformer patterns, custom transformers	Reference
Evidently AI (evidently.ai)	Library	Open-source drift monitoring, production-ready	2 hrs
“ML Engineering” — Andriy Burkov	Book	Practical deployment patterns, no fluff	6 hrs

---

Return to README.md · Previous: 04_feature_engineering.md · Next: 06_phase_projects.md