Logistic Regression

What this module covers:

import numpy as np
import pandas as pd
import matplotlib
matplotlib.use('Agg')   # non-interactive backend for scripts

# Demonstrate why linear regression fails for classification
np.random.seed(42)
n = 200

# Generate binary classification data: late vs on-time deliveries
distance_ontime = np.random.normal(3, 1, n//2)
distance_late   = np.random.normal(6, 1, n//2)

X = np.concatenate([distance_ontime, distance_late]).reshape(-1, 1)
y = np.concatenate([np.zeros(n//2), np.ones(n//2)])   # 0=on-time, 1=late

from sklearn.linear_model import LinearRegression, LogisticRegression
from sklearn.preprocessing import StandardScaler

sc = StandardScaler()
X_sc = sc.fit_transform(X)

# Linear regression on binary labels
lin_reg = LinearRegression()
lin_reg.fit(X_sc, y)
y_pred_lin = lin_reg.predict(X_sc)

print("Linear regression predictions (should be in [0,1] but are not):")
print(f"  Min prediction: {y_pred_lin.min():.3f}")
print(f"  Max prediction: {y_pred_lin.max():.3f}")
print(f"  Predictions < 0: {(y_pred_lin < 0).sum()}")
print(f"  Predictions > 1: {(y_pred_lin > 1).sum()}")

# Logistic regression
log_reg = LogisticRegression()
log_reg.fit(X_sc, y)
y_pred_proba = log_reg.predict_proba(X_sc)[:, 1]

print("
Logistic regression predictions (always in [0,1]):")
print(f"  Min prediction: {y_pred_proba.min():.3f}")
print(f"  Max prediction: {y_pred_proba.max():.3f}")
print(f"  Predictions < 0: {(y_pred_proba < 0).sum()}")
print(f"  Predictions > 1: {(y_pred_proba > 1).sum()}")

import numpy as np

# ── The sigmoid function ──────────────────────────────────────────────
def sigmoid(z: np.ndarray) -> np.ndarray:
    """σ(z) = 1 / (1 + e^(-z)) — numerically stable implementation."""
    # Clip to prevent overflow in exp for very large negative values
    return 1 / (1 + np.exp(-np.clip(z, -500, 500)))

# Properties of sigmoid
z_values = np.array([-10, -5, -2, -1, 0, 1, 2, 5, 10])
print("Sigmoid values:")
for z, s in zip(z_values, sigmoid(z_values)):
    print(f"  σ({z:+3d}) = {s:.6f}")

# Sigmoid derivative: σ'(z) = σ(z) * (1 - σ(z))
# Maximum at z=0: σ'(0) = 0.5 * 0.5 = 0.25
# This is why sigmoid causes vanishing gradients in deep networks

def sigmoid_derivative(z):
    s = sigmoid(z)
    return s * (1 - s)

print(f"
Sigmoid derivative at z=0: {sigmoid_derivative(0):.4f}")  # 0.25 (max)
print(f"Sigmoid derivative at z=5: {sigmoid_derivative(5):.6f}")   # near 0

# ── Logistic regression prediction ───────────────────────────────────
def logistic_predict(X: np.ndarray, w: np.ndarray, b: float) -> np.ndarray:
    """
    Forward pass: compute probability for each sample.
    X: (n_samples, n_features)
    w: (n_features,)
    b: scalar
    Returns: (n_samples,) — probability of class 1
    """
    z = X @ w + b         # linear combination: shape (n_samples,)
    return sigmoid(z)     # squash to [0,1]

# Test with one sample: distance=6km, traffic=8, prep=20min
# High values → should predict high probability of being late
X_sample = np.array([[6.0, 8.0, 20.0]])   # shape (1, 3)
w_random  = np.array([0.5, 0.3, 0.2])     # random weights
b_random  = -5.0

p_late = logistic_predict(X_sample, w_random, b_random)
print(f"
P(late | distance=6, traffic=8, prep=20) = {p_late[0]:.4f}")
print(f"Prediction: {'LATE' if p_late[0] >= 0.5 else 'ON-TIME'}")

import numpy as np

def binary_cross_entropy(y_true: np.ndarray, y_pred: np.ndarray,
                          eps: float = 1e-10) -> float:
    """
    Binary cross-entropy loss.
    eps prevents log(0) = -inf.
    """
    y_pred = np.clip(y_pred, eps, 1 - eps)   # prevent log(0)
    return -np.mean(
        y_true * np.log(y_pred) + (1 - y_true) * np.log(1 - y_pred)
    )

# Loss for different prediction / label combinations
scenarios = [
    (1, 0.99,  'Truly late,   predicts 0.99 (confident & correct)'),
    (1, 0.50,  'Truly late,   predicts 0.50 (uncertain)'),
    (1, 0.01,  'Truly late,   predicts 0.01 (confident & WRONG)'),
    (0, 0.01,  'Truly on-time, predicts 0.01 (confident & correct)'),
    (0, 0.50,  'Truly on-time, predicts 0.50 (uncertain)'),
    (0, 0.99,  'Truly on-time, predicts 0.99 (confident & WRONG)'),
]

print("Cross-entropy loss per scenario:")
for y, p, desc in scenarios:
    loss = binary_cross_entropy(np.array([y]), np.array([p]))
    bar  = '█' * int(loss * 5)
    print(f"  {desc}")
    print(f"    Loss = {loss:.4f}  {bar}")
    print()

# Key insight: loss is SYMMETRIC
# L(y=1, p=0.01) ≈ L(y=0, p=0.99) ≈ 4.6
# Both confident wrong predictions are equally penalised

import numpy as np

np.random.seed(42)

# ── Generate DoorDash dataset ───────────────────────────────────────────
n = 2000
distance = np.abs(np.random.normal(4.0, 2.0, n)).clip(0.5, 15)
traffic  = np.random.randint(1, 11, n).astype(float)
prep     = np.abs(np.random.normal(15, 5, n)).clip(5, 35)
delivery = 8.6 + 7.3*distance + 0.8*prep + 1.5*traffic + np.random.normal(0, 4, n)
y = (delivery > 45).astype(float)   # 1 = late, 0 = on-time

X = np.column_stack([distance, traffic, prep])

# Scale features (critical for gradient descent convergence)
from sklearn.preprocessing import StandardScaler
from sklearn.model_selection import train_test_split

X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)
sc = StandardScaler()
X_train_sc = sc.fit_transform(X_train)
X_test_sc  = sc.transform(X_test)

# ── Logistic regression from scratch ─────────────────────────────────
class LogisticRegressionScratch:
    """
    Binary logistic regression trained with mini-batch gradient descent.
    All the math visible, nothing hidden.
    """

    def __init__(self, lr: float = 0.1, n_epochs: int = 200,
                 batch_size: int = 64, l2: float = 0.01):
        self.lr          = lr
        self.n_epochs    = n_epochs
        self.batch_size  = batch_size
        self.l2          = l2   # L2 regularisation strength
        self.w           = None
        self.b           = 0.0
        self.loss_history = []

    @staticmethod
    def sigmoid(z): return 1 / (1 + np.exp(-np.clip(z, -500, 500)))

    def fit(self, X: np.ndarray, y: np.ndarray):
        n, d = X.shape
        self.w = np.zeros(d)   # initialise weights at 0
        rng    = np.random.default_rng(42)

        for epoch in range(self.n_epochs):
            # Shuffle training data each epoch
            idx = rng.permutation(n)
            X_s, y_s = X[idx], y[idx]
            epoch_loss = 0.0

            for start in range(0, n, self.batch_size):
                Xb = X_s[start:start + self.batch_size]
                yb = y_s[start:start + self.batch_size]
                nb = len(Xb)

                # ── Forward pass ──────────────────────────────────
                z = Xb @ self.w + self.b        # linear: (nb,)
                p = self.sigmoid(z)              # probability: (nb,)

                # ── Loss (with L2 regularisation) ─────────────────
                eps   = 1e-10
                bce   = -np.mean(yb*np.log(p+eps) + (1-yb)*np.log(1-p+eps))
                reg   = 0.5 * self.l2 * np.dot(self.w, self.w)
                loss  = bce + reg
                epoch_loss += loss

                # ── Backward pass (gradients) ─────────────────────
                error  = p - yb                        # (nb,) — prediction error
                dw     = (Xb.T @ error) / nb + self.l2 * self.w  # (d,)
                db     = error.mean()                  # scalar

                # ── Weight update ──────────────────────────────────
                self.w -= self.lr * dw
                self.b -= self.lr * db

            self.loss_history.append(epoch_loss / (n // self.batch_size))

            if epoch % 40 == 0:
                print(f"  Epoch {epoch:3d}: loss = {self.loss_history[-1]:.4f}")

    def predict_proba(self, X: np.ndarray) -> np.ndarray:
        return self.sigmoid(X @ self.w + self.b)

    def predict(self, X: np.ndarray, threshold: float = 0.5) -> np.ndarray:
        return (self.predict_proba(X) >= threshold).astype(int)

    def score(self, X: np.ndarray, y: np.ndarray) -> float:
        return (self.predict(X) == y).mean()

# ── Train ──────────────────────────────────────────────────────────────
model_scratch = LogisticRegressionScratch(lr=0.1, n_epochs=200, l2=0.01)
model_scratch.fit(X_train_sc, y_train)

acc_train = model_scratch.score(X_train_sc, y_train)
acc_test  = model_scratch.score(X_test_sc, y_test)
print(f"
From-scratch model:")
print(f"  Train accuracy: {acc_train:.4f}")
print(f"  Test accuracy:  {acc_test:.4f}")

# Verify against sklearn
from sklearn.linear_model import LogisticRegression
sk_model = LogisticRegression(C=1/0.01, max_iter=500, random_state=42)
sk_model.fit(X_train_sc, y_train)
print(f"
sklearn LogisticRegression test accuracy: {sk_model.score(X_test_sc, y_test):.4f}")
# Should match closely ↑

from sklearn.linear_model import LogisticRegression
from sklearn.pipeline import Pipeline
from sklearn.preprocessing import StandardScaler
from sklearn.model_selection import cross_val_score
import numpy as np

# ── The C parameter — inverse of regularisation strength ──────────────
# C = 1/lambda  where lambda is the regularisation strength
# Large C  → weak regularisation → model can fit training data more freely
# Small C  → strong regularisation → simpler model, less overfitting
# Default C=1.0 is a reasonable starting point
# ALWAYS tune C — it is the most important hyperparameter

for C in [0.001, 0.01, 0.1, 1.0, 10.0, 100.0]:
    model = LogisticRegression(C=C, max_iter=1000, random_state=42)
    scores = cross_val_score(model, X_train_sc, y_train, cv=5, scoring='roc_auc')
    print(f"  C={C:<8}: ROC-AUC = {scores.mean():.4f} ± {scores.std():.4f}")

# ── penalty — L1 vs L2 vs ElasticNet ──────────────────────────────────
# penalty='l2' (default): weight² penalty — drives weights toward 0, keeps all features
# penalty='l1':           |weight| penalty — drives some weights to exactly 0 (feature selection)
# penalty='elasticnet':   mix of L1 and L2 — l1_ratio controls the mix

penalties = [
    ('l2',         dict(C=1.0, solver='lbfgs')),
    ('l1',         dict(C=1.0, solver='liblinear')),
    ('elasticnet', dict(C=1.0, solver='saga', l1_ratio=0.5, max_iter=2000)),
    ('none',       dict(solver='lbfgs')),
]

print("
Regularisation comparison (5-fold CV ROC-AUC):")
for penalty, kwargs in penalties:
    model = LogisticRegression(penalty=penalty, random_state=42, **kwargs)
    scores = cross_val_score(model, X_train_sc, y_train, cv=5, scoring='roc_auc')
    print(f"  penalty='{penalty}'   : {scores.mean():.4f} ± {scores.std():.4f}")

# ── solver — which optimisation algorithm ─────────────────────────────
# 'lbfgs':     default, L-BFGS quasi-Newton. Fast for small-medium data.
# 'liblinear': fast for small datasets, supports L1 and L2
# 'saga':      faster for large datasets, supports L1, L2, ElasticNet
# 'sag':       variant of SAGA, only L2
# 'newton-cg': Newton's method, only L2

# ── class_weight — handle imbalanced datasets ─────────────────────────
# If 90% of orders are on-time and 10% are late, a model that always
# predicts on-time gets 90% accuracy but catches 0 late orders.
# class_weight='balanced' weights the loss by inverse class frequency.

model_balanced = LogisticRegression(
    C=1.0,
    class_weight='balanced',   # auto-weight: rare class gets higher penalty for errors
    max_iter=1000,
    random_state=42,
)
model_balanced.fit(X_train_sc, y_train)

from sklearn.metrics import classification_report
y_pred_balanced = model_balanced.predict(X_test_sc)
print("
With class_weight='balanced':")
print(classification_report(y_test, y_pred_balanced, target_names=['on-time','late']))

# ── Full pipeline — the production pattern ─────────────────────────────
pipe = Pipeline([
    ('scaler', StandardScaler()),
    ('model',  LogisticRegression(
        C=1.0,
        penalty='l2',
        solver='lbfgs',
        max_iter=1000,
        class_weight='balanced',
        random_state=42,
    )),
])
pipe.fit(X_train, y_train)
print(f"
Pipeline test accuracy: {pipe.score(X_test, y_test):.4f}")

import numpy as np
import pandas as pd
from sklearn.linear_model import LogisticRegression
from sklearn.preprocessing import StandardScaler

# Train on full feature set including engineered features
feature_names = ['distance_km','traffic_score','restaurant_prep']
sc = StandardScaler()
X_train_sc = sc.fit_transform(X_train)
X_test_sc  = sc.transform(X_test)

model = LogisticRegression(C=1.0, max_iter=1000, random_state=42)
model.fit(X_train_sc, y_train)

# ── Coefficient interpretation ────────────────────────────────────────
# Coefficients are in terms of STANDARDISED features
# A coefficient of 1.5 means: increasing this feature by 1 std
# increases the log-odds of being late by 1.5
# log-odds = log(p / (1-p))

coef_df = pd.DataFrame({
    'feature':     feature_names,
    'coefficient': model.coef_[0],
    'odds_ratio':  np.exp(model.coef_[0]),  # e^coef = how odds multiply per std
}).sort_values('coefficient', ascending=False)

print("Logistic regression coefficients (standardised features):")
print(coef_df.to_string(index=False))
print(f"
Intercept (bias): {model.intercept_[0]:.4f}")

print("
Interpretation (per 1 standard deviation increase):")
for _, row in coef_df.iterrows():
    direction = 'INCREASES' if row['coefficient'] > 0 else 'DECREASES'
    print(f"  {row['feature']:<20}: {direction} P(late) | "
          f"coef={row['coefficient']:.3f} | odds_ratio={row['odds_ratio']:.3f}")

# ── Decision boundary in 2D (distance vs traffic) ─────────────────────
# Boundary: w₁x₁ + w₂x₂ + b = 0  →  x₂ = -(w₁x₁ + b) / w₂
# For a 2-feature model:
model_2d = LogisticRegression(C=1.0, random_state=42)
X_2d_train = X_train_sc[:, :2]   # only distance and traffic (standardised)
model_2d.fit(X_2d_train, y_train)

w1, w2 = model_2d.coef_[0]
b      = model_2d.intercept_[0]

# x_traffic = -(w1*x_distance + b) / w2
print(f"
Decision boundary equation:")
print(f"  {w2:.3f} × traffic + {w1:.3f} × distance + {b:.3f} = 0")
print(f"  traffic = ({-w1:.3f} × distance + {-b:.3f}) / {w2:.3f}")

# At distance = 0 (standardised): traffic threshold
traffic_at_mean_dist = (- w1 * 0 - b) / w2
print(f"  At mean distance: classify as late if traffic > {traffic_at_mean_dist:.2f} std")

from sklearn.metrics import (
    classification_report, confusion_matrix,
    roc_auc_score, average_precision_score,
    precision_recall_curve, roc_curve,
)
import numpy as np

# Train final model
model = LogisticRegression(C=1.0, class_weight='balanced',
                            max_iter=1000, random_state=42)
model.fit(X_train_sc, y_train)

y_pred       = model.predict(X_test_sc)
y_proba      = model.predict_proba(X_test_sc)[:, 1]

# ── Classification report ──────────────────────────────────────────────
print("Classification report:")
print(classification_report(y_test, y_pred, target_names=['on-time','late']))

# ── Confusion matrix ───────────────────────────────────────────────────
cm = confusion_matrix(y_test, y_pred)
tn, fp, fn, tp = cm.ravel()
print(f"Confusion matrix:")
print(f"  TN={tn}  FP={fp}")
print(f"  FN={fn}  TP={tp}")
print(f"
  Precision (of predicted-late, how many truly late): {tp/(tp+fp):.3f}")
print(f"  Recall    (of truly-late, how many we caught):       {tp/(tp+fn):.3f}")

# ── ROC-AUC ───────────────────────────────────────────────────────────
roc_auc = roc_auc_score(y_test, y_proba)
avg_prec = average_precision_score(y_test, y_proba)
print(f"
  ROC-AUC:              {roc_auc:.4f}")
print(f"  Avg Precision (AP):   {avg_prec:.4f}")
print(f"  Accuracy:             {(y_pred == y_test).mean():.4f}")

# ── Threshold tuning — 0.5 is not always optimal ──────────────────────
# The business question determines the right threshold.
# At DoorDash: missing a late delivery (FN) costs customer experience.
# A lower threshold catches more late deliveries but flags more on-time ones.

print("
Threshold analysis:")
print(f"{'Threshold':<12} {'Precision':<12} {'Recall':<10} {'F1':<10} {'Flagged %'}")
print("─" * 58)
for threshold in [0.2, 0.3, 0.4, 0.5, 0.6, 0.7]:
    y_pred_t  = (y_proba >= threshold).astype(int)
    tp_t      = ((y_pred_t == 1) & (y_test == 1)).sum()
    fp_t      = ((y_pred_t == 1) & (y_test == 0)).sum()
    fn_t      = ((y_pred_t == 0) & (y_test == 1)).sum()
    prec_t    = tp_t / max(tp_t + fp_t, 1)
    rec_t     = tp_t / max(tp_t + fn_t, 1)
    f1_t      = 2 * prec_t * rec_t / max(prec_t + rec_t, 1e-10)
    flagged_t = y_pred_t.mean() * 100
    print(f"{threshold:<12.1f} {prec_t:<12.3f} {rec_t:<10.3f} {f1_t:<10.3f} {flagged_t:.1f}%")

# For DoorDash: a threshold of 0.3 catches more late deliveries
# but flags more on-time ones for proactive ETA warnings — good tradeoff

import numpy as np
import pandas as pd
from sklearn.linear_model import LogisticRegression
from sklearn.preprocessing import StandardScaler

# Generate a high-dimensional problem: 20 features, many irrelevant
np.random.seed(42)
n_train = 500

# 3 truly useful features + 17 noise features
X_useful = X_train[:n_train, :]                            # (n, 3) useful
X_noise  = np.random.randn(n_train, 17)                   # (n, 17) pure noise
X_high   = np.column_stack([X_useful, X_noise])
y_high   = y_train[:n_train]

X_useful_te = X_test[:100, :]
X_noise_te  = np.random.randn(100, 17)
X_high_te   = np.column_stack([X_useful_te, X_noise_te])
y_high_te   = y_test[:100]

sc20 = StandardScaler()
X_high_sc    = sc20.fit_transform(X_high)
X_high_te_sc = sc20.transform(X_high_te)

feature_names_20 = (
    ['distance_km','traffic_score','restaurant_prep'] +
    [f'noise_{i:02d}' for i in range(17)]
)

print("Effect of regularisation on 20-feature problem (3 useful, 17 noise):")
print(f"{'Model':<35} {'Train acc':<12} {'Test acc':<12} {'Non-zero coefs'}")
print("─" * 72)

for name, model in [
    ('No regularisation (penalty=None)', LogisticRegression(penalty=None, solver='lbfgs', max_iter=2000)),
    ('L2 regularisation C=1.0',          LogisticRegression(penalty='l2', C=1.0, solver='lbfgs', max_iter=2000)),
    ('L2 regularisation C=0.1',          LogisticRegression(penalty='l2', C=0.1, solver='lbfgs', max_iter=2000)),
    ('L1 regularisation C=1.0',          LogisticRegression(penalty='l1', C=1.0, solver='liblinear', max_iter=2000)),
    ('L1 regularisation C=0.1',          LogisticRegression(penalty='l1', C=0.1, solver='liblinear', max_iter=2000)),
]:
    model.fit(X_high_sc, y_high)
    tr_acc = model.score(X_high_sc, y_high)
    te_acc = model.score(X_high_te_sc, y_high_te)
    n_nz   = (np.abs(model.coef_[0]) > 1e-6).sum()
    print(f"{name:<35} {tr_acc:<12.4f} {te_acc:<12.4f} {n_nz}/20")

# L1 observations:
# - Drives noise feature coefficients to exactly 0 (true sparse solution)
# - Keeps only the 3 useful features (if C is small enough)
# - Acts as automatic feature selection

# L2 observations:
# - Shrinks all coefficients toward 0 but never exactly 0
# - All 20 features stay in the model with small weights
# - Better when all features are somewhat relevant

# Check which features L1 kept
l1_model = LogisticRegression(penalty='l1', C=0.1, solver='liblinear', max_iter=2000)
l1_model.fit(X_high_sc, y_high)
kept = [name for name, coef in zip(feature_names_20, l1_model.coef_[0]) if abs(coef) > 1e-6]
print(f"
Features kept by L1 (C=0.1): {kept}")

from sklearn.linear_model import LogisticRegression
from sklearn.preprocessing import StandardScaler, LabelEncoder
from sklearn.metrics import classification_report
import numpy as np
import pandas as pd

np.random.seed(42)
n = 3000

# Multi-class problem: predict delivery speed category
# Classes: 'express' (<25 min), 'normal' (25-45 min), 'delayed' (>45 min)
distance = np.abs(np.random.normal(4.0, 2.0, n)).clip(0.5, 15)
traffic  = np.random.randint(1, 11, n).astype(float)
prep     = np.abs(np.random.normal(15, 5, n)).clip(5, 35)
delivery = 8.6 + 7.3*distance + 0.8*prep + 1.5*traffic + np.random.normal(0, 4, n)

y_multi = pd.cut(
    delivery,
    bins=[-np.inf, 25, 45, np.inf],
    labels=['express','normal','delayed'],
).astype(str)

le = LabelEncoder()
y_enc = le.fit_transform(y_multi)
print(f"Classes: {le.classes_}")
print(f"Class distribution: {dict(zip(le.classes_, np.bincount(y_enc)))}")

X_multi = np.column_stack([distance, traffic, prep])
from sklearn.model_selection import train_test_split
Xm_tr, Xm_te, ym_tr, ym_te = train_test_split(X_multi, y_enc,
                                                test_size=0.2, random_state=42)
sc = StandardScaler()
Xm_tr_sc = sc.fit_transform(Xm_tr)
Xm_te_sc = sc.transform(Xm_te)

# ── OvR ────────────────────────────────────────────────────────────────
model_ovr = LogisticRegression(
    multi_class='ovr',
    C=1.0, solver='lbfgs', max_iter=1000, random_state=42
)
model_ovr.fit(Xm_tr_sc, ym_tr)

# ── Multinomial (Softmax) ──────────────────────────────────────────────
model_softmax = LogisticRegression(
    multi_class='multinomial',
    C=1.0, solver='lbfgs', max_iter=1000, random_state=42
)
model_softmax.fit(Xm_tr_sc, ym_tr)

print(f"
OvR accuracy:       {model_ovr.score(Xm_te_sc, ym_te):.4f}")
print(f"Softmax accuracy:   {model_softmax.score(Xm_te_sc, ym_te):.4f}")

print("
Softmax probabilities for one sample:")
sample = Xm_te_sc[:1]
proba  = model_softmax.predict_proba(sample)[0]
for cls, p in zip(le.classes_, proba):
    bar = '█' * int(p * 30)
    print(f"  {cls:<10}: {bar} {p:.4f}")
print(f"  Sum: {proba.sum():.6f}  ← always exactly 1.0")

print("
Classification report (Softmax):")
print(classification_report(
    ym_te, model_softmax.predict(Xm_te_sc),
    target_names=le.classes_,
))

import numpy as np
import pandas as pd
from sklearn.pipeline import Pipeline
from sklearn.compose import ColumnTransformer
from sklearn.preprocessing import StandardScaler, OneHotEncoder
from sklearn.impute import SimpleImputer
from sklearn.linear_model import LogisticRegression
from sklearn.model_selection import StratifiedKFold, cross_validate
from sklearn.metrics import (roc_auc_score, average_precision_score,
                              classification_report)
import joblib

np.random.seed(42)
n = 8000
restaurants = ['Pizza Hut','Biryani Blues',"McDonald's","Haldiram's",
               'Dominos','KFC','Subway','Burger King']
cities = ['Seattle','New York','Denver','Austin','Boston','Chicago']
slots  = ['breakfast','lunch','evening','dinner']

distance = np.abs(np.random.normal(4.0, 2.0, n)).clip(0.5, 15)
traffic  = np.random.randint(1, 11, n).astype(float)
prep     = np.abs(np.random.normal(15, 5, n)).clip(5, 35)
value    = np.abs(np.random.normal(350, 150, n)).clip(50, 1200)
delivery = (8.6 + 7.3*distance + 0.8*prep + 1.5*traffic
            + np.random.normal(0, 4, n)).clip(10, 120)

df_prod = pd.DataFrame({
    'restaurant':     np.random.choice(restaurants, n),
    'city':           np.random.choice(cities, n),
    'time_slot':      np.random.choice(slots, n),
    'distance_km':    distance,
    'traffic_score':  traffic,
    'restaurant_prep':prep,
    'order_value':    value,
})
y_prod = (delivery > 45).astype(int)

# ── Feature engineering ────────────────────────────────────────────────
df_prod['log_distance']  = np.log1p(df_prod['distance_km'])
df_prod['dist_x_traffic']= df_prod['distance_km'] * df_prod['traffic_score']
df_prod['log_value']     = np.log1p(df_prod['order_value'])

NUM_FEATURES = ['log_distance','traffic_score','restaurant_prep',
                'dist_x_traffic','log_value','distance_km']
CAT_FEATURES = ['city','time_slot']

# ── Column transformer ─────────────────────────────────────────────────
preprocessor = ColumnTransformer([
    ('num', Pipeline([
        ('imputer', SimpleImputer(strategy='median')),
        ('scaler',  StandardScaler()),
    ]), NUM_FEATURES),
    ('cat', Pipeline([
        ('imputer', SimpleImputer(strategy='most_frequent')),
        ('onehot',  OneHotEncoder(handle_unknown='ignore', sparse_output=False, drop='first')),
    ]), CAT_FEATURES),
])

# ── Final pipeline ─────────────────────────────────────────────────────
pipeline = Pipeline([
    ('preprocessor', preprocessor),
    ('classifier',   LogisticRegression(
        C=0.5,
        penalty='l2',
        class_weight='balanced',
        solver='lbfgs',
        max_iter=2000,
        random_state=42,
    )),
])

# ── Stratified cross-validation ────────────────────────────────────────
cv = StratifiedKFold(n_splits=5, shuffle=True, random_state=42)
cv_results = cross_validate(
    pipeline, df_prod, y_prod,
    cv=cv,
    scoring=['roc_auc','average_precision','accuracy'],
    return_train_score=True,
)

print("5-Fold Cross-Validation Results:")
for metric in ['roc_auc','average_precision','accuracy']:
    train_mean = cv_results[f'train_{metric}'].mean()
    val_mean   = cv_results[f'test_{metric}'].mean()
    val_std    = cv_results[f'test_{metric}'].std()
    print(f"  {metric:<22}: train={train_mean:.4f}  val={val_mean:.4f} ± {val_std:.4f}")

# ── Train final model on full data ─────────────────────────────────────
pipeline.fit(df_prod, y_prod)

# ── Find optimal threshold on held-out validation set ─────────────────
from sklearn.model_selection import train_test_split
df_tr, df_val, y_tr, y_val = train_test_split(df_prod, y_prod,
                                               test_size=0.2, stratify=y_prod, random_state=99)
pipeline.fit(df_tr, y_tr)
val_proba = pipeline.predict_proba(df_val)[:, 1]

# Choose threshold that maximises F1 for the late class
from sklearn.metrics import f1_score
best_threshold, best_f1 = 0.5, 0.0
for t in np.arange(0.2, 0.8, 0.02):
    f1 = f1_score(y_val, (val_proba >= t).astype(int), pos_label=1)
    if f1 > best_f1:
        best_f1, best_threshold = f1, t

print(f"
Optimal threshold: {best_threshold:.2f} (F1={best_f1:.4f})")

# ── Save the pipeline ──────────────────────────────────────────────────
joblib.dump({
    'pipeline':  pipeline,
    'threshold': best_threshold,
    'features':  NUM_FEATURES + CAT_FEATURES,
    'version':   'v1.0',
}, '/tmp/late_delivery_model.pkl')
print("Model saved to /tmp/late_delivery_model.pkl")

# ── Load and score a new order ─────────────────────────────────────────
saved = joblib.load('/tmp/late_delivery_model.pkl')
new_order = pd.DataFrame([{
    'restaurant': 'Pizza Hut', 'city': 'Seattle', 'time_slot': 'dinner',
    'distance_km': 7.5, 'traffic_score': 9, 'restaurant_prep': 22,
    'order_value': 480,
}])
new_order['log_distance']   = np.log1p(new_order['distance_km'])
new_order['dist_x_traffic'] = new_order['distance_km'] * new_order['traffic_score']
new_order['log_value']      = np.log1p(new_order['order_value'])

p_late = saved['pipeline'].predict_proba(new_order)[0, 1]
pred   = 'LATE' if p_late >= saved['threshold'] else 'ON-TIME'
print(f"
New order prediction: {pred} (P(late)={p_late:.3f})")

Swapping in the sigmoid is not a decorative final step — the entire loss function had to change to make the combination well-behaved. Minimising squared error through a sigmoid produces a non-convex, multi-modal loss surface that gradient descent can get stuck in. Logistic regression instead pairs the sigmoid with cross-entropy loss specifically because that combination is provably convex, guaranteeing gradient descent finds the global minimum. The sigmoid is not bolted onto linear regression's loss; the loss itself had to change for the model to work reliably.

Coefficients live in log-odds space, not probability space. A coefficient of 1.5 does not mean "adds 1.5 to the probability" — it means a one-unit increase in that feature multiplies the odds of the positive class by e^1.5, about 4.5x. Odds and probability are related but different things (odds equals p divided by one minus p), and turning a coefficient into an actual probability change requires running the whole linear combination back through the sigmoid. The effect of one feature on probability also depends on the values of every other feature — it is not a constant additive effect the way it would be in ordinary linear regression.

It trades one set of assumptions for a different, still-real set. Instead of assuming a linear relationship between features and the outcome directly, it assumes a linear relationship between features and the log-odds of the outcome — meaning the decision boundary it can draw is still a straight line, or hyperplane, in feature space, just wrapped in a sigmoid. It still assumes independent observations, needs the absence of severe multicollinearity for stable coefficients, and wants a reasonably large sample per feature to avoid unstable estimates. Fewer assumptions is not the same thing as no assumptions.

They solve different problems. L2 shrinks every coefficient toward zero but essentially never sets one to exactly zero — the right choice when most features probably carry some real signal and the goal is just to tame their magnitude. L1 can drive coefficients to exactly zero, performing automatic feature selection — better when many features are suspected to be pure noise and a sparse, interpretable model is preferred. Picking the wrong one either leaves noisy features in the model at full strength (L2 on a noisy feature set) or arbitrarily discards useful correlated features (L1 tends to keep one from a correlated group and zero out the rest).

0.5 is only correct when a false positive and a false negative cost exactly the same, which is rare in practice. A fraud model where missing real fraud is far more costly than a false alarm should use a lower threshold, flagging more transactions to catch more fraud while accepting more false alarms. A screening test with expensive or invasive follow-ups might want a higher threshold instead. The threshold is a business decision made by weighing the cost of each error type, not a property the model hands you — that is exactly why ROC and precision-recall curves exist: to let that point be chosen deliberately.

MSE paired with a sigmoid produces a non-convex loss surface — because the sigmoid itself is non-linear, squaring the error re-introduces multiple local minima that gradient descent can get stuck in, with no guarantee of finding the best fit. Cross-entropy, paired with sigmoid, produces a convex loss surface, so gradient descent is guaranteed to converge to the global minimum. There is also a modelling reason: cross-entropy penalises confident wrong predictions far more harshly (predicting 0.99 when the true label is 0 costs about 4.6, versus MSE's roughly 0.98), which teaches the model to only be confident when it is actually likely to be correct — exactly the behaviour you want from a probability estimator.

A one standard deviation increase in income multiplies the odds of the positive outcome by e^0.7, roughly 2.0x, holding every other feature constant. I would not say it increases the probability by 70%, or by any fixed amount — the actual probability shift depends on where you start on the sigmoid curve: near the middle, around probability 0.5, a coefficient like this shifts probability substantially; near the tails, close to 0 or 1, the same coefficient barely moves probability at all because the sigmoid saturates there. I would also check whether income was actually standardised, since the coefficient is only comparable to the others if every feature shares the same scale.

Left at defaults, the model can hit 98% accuracy by predicting "not fraud" for every single transaction, since that already matches the base rate — the loss function has little incentive to fit the rare class well. I would stop trusting accuracy entirely and evaluate with precision, recall, and ROC-AUC or PR-AUC instead. For the training itself, I would use class_weight='balanced' so errors on the minority class are weighted more heavily in the loss, or oversample or undersample as an alternative. Just as important, I would not leave the decision threshold at 0.5 — I would pick a threshold on a validation set that reflects the real cost of missing fraud against the cost of a false alarm, since that is usually a bigger lever than the model architecture itself.

Probability is bounded between 0 and 1 and is what most people intuitively reason about. Odds is probability divided by one minus probability, ranging from 0 up to infinity — unbounded on the upper end. Logistic regression's linear component can output any real number, and the log of the odds is the one transformation of probability that also ranges over every real number, which is exactly why the model is linear in log-odds space rather than linear in probability space directly. Coefficients are reported as effects on log-odds, or on odds after exponentiating, because that is the space where the model's math is genuinely linear; reporting them as probability effects would be misleading, since that relationship is neither linear nor constant.

When the true decision boundary is clearly non-linear even in log-odds space — say the relationship between a feature and the outcome is U-shaped — plain logistic regression will underfit unless polynomial or interaction terms are hand-engineered in. I would also avoid it when I need complex feature interactions captured automatically, since tree-based models like XGBoost pick those up natively while logistic regression requires them to be specified by hand. And if calibrated probabilities are not actually needed, gradient boosting methods usually beat logistic regression on raw predictive accuracy for complex tabular data out of the box. I would still often start with logistic regression anyway, though, as a fast, interpretable baseline before reaching for something heavier.