S086 · Logistic Regression — Sigmoid, Cross-Entropy, from Scratch

# logreg_scratch.py — logistic regression, ~130 lines runnable.
import numpy as np
from sklearn.datasets import load_breast_cancer
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler
from sklearn.linear_model import LogisticRegression
from sklearn.metrics import (
    accuracy_score, precision_recall_fscore_support, roc_auc_score,
    precision_recall_curve, average_precision_score, confusion_matrix
)
 
np.random.seed(42)
 
# --- 1. Data + split ---
data = load_breast_cancer()
X, y = data.data, data.target        # y=1 → benign, y=0 → malignant
print(f"X shape: {X.shape}, class balance: y=1 → {y.mean():.2%}")
 
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, stratify=y, random_state=42)
 
scaler = StandardScaler().fit(X_train)
X_train_s = scaler.transform(X_train)
X_test_s  = scaler.transform(X_test)
 
def add_bias(X):
    return np.hstack([np.ones((X.shape[0], 1)), X])
 
Xtr = add_bias(X_train_s)
Xte = add_bias(X_test_s)
 
# --- 2. Sigmoid + cross-entropy ---
def sigmoid(z):
    """Numerically stable sigmoid. Avoids overflow for very large |z|."""
    return np.where(z >= 0, 1 / (1 + np.exp(-z)), np.exp(z) / (1 + np.exp(z)))
 
def bce_loss(y, p, eps=1e-15):
    """Binary cross-entropy. Clip p to avoid log(0)."""
    p = np.clip(p, eps, 1 - eps)
    return -np.mean(y * np.log(p) + (1 - y) * np.log(1 - p))
 
# --- 3. Fit via gradient descent ---
def fit_logreg(X, y, lr=0.1, epochs=2000, verbose=False):
    n, d = X.shape
    theta = np.zeros(d)
    for epoch in range(epochs):
        z = X @ theta
        p = sigmoid(z)
        grad = (1 / n) * (X.T @ (p - y))           # ← identical shape to linreg gradient
        theta = theta - lr * grad
        if verbose and epoch % 200 == 0:
            print(f"  epoch {epoch:4d}: loss={bce_loss(y, p):.4f}, acc={(p > 0.5).astype(int).mean():.3f}")
    return theta
 
theta = fit_logreg(Xtr, y_train, lr=0.1, epochs=2000, verbose=True)
 
# --- 4. Predict + evaluate ---
p_test = sigmoid(Xte @ theta)
yhat_default = (p_test > 0.5).astype(int)
 
acc = accuracy_score(y_test, yhat_default)
auc = roc_auc_score(y_test, p_test)
pr_auc = average_precision_score(y_test, p_test)
print(f"\n[scratch @ threshold 0.5]  acc={acc:.4f}  ROC-AUC={auc:.4f}  PR-AUC={pr_auc:.4f}")
 
# --- 5. sklearn sanity check ---
sk = LogisticRegression(penalty=None, max_iter=5000).fit(X_train_s, y_train)
p_sk = sk.predict_proba(X_test_s)[:, 1]
print(f"[sklearn @ threshold 0.5]  acc={accuracy_score(y_test, sk.predict(X_test_s)):.4f}  "
      f"ROC-AUC={roc_auc_score(y_test, p_sk):.4f}")
 
# --- 6. Threshold tuning via PR curve ---
prec, rec, thr = precision_recall_curve(y_test, p_test)
f1 = 2 * prec * rec / (prec + rec + 1e-12)
best_idx = np.argmax(f1[:-1])         # skip last point (no threshold)
best_thr = thr[best_idx]
print(f"\nBest F1 = {f1[best_idx]:.4f} at threshold = {best_thr:.4f}")
print(f"  → precision = {prec[best_idx]:.4f}, recall = {rec[best_idx]:.4f}")
 
yhat_tuned = (p_test > best_thr).astype(int)
print(f"[scratch @ tuned threshold] confusion matrix:\n{confusion_matrix(y_test, yhat_tuned)}")
 
# --- 7. Calibration check — reliability plot data ---
def reliability(y, p, n_bins=10):
    """Bin predictions by confidence, check if the fraction of positives matches."""
    bins = np.linspace(0, 1, n_bins + 1)
    idx = np.digitize(p, bins) - 1
    idx = np.clip(idx, 0, n_bins - 1)
    means_pred, means_true, counts = [], [], []
    for b in range(n_bins):
        mask = idx == b
        if mask.sum() > 0:
            means_pred.append(p[mask].mean())
            means_true.append(y[mask].mean())
            counts.append(int(mask.sum()))
    return means_pred, means_true, counts
 
mp, mt, ct = reliability(y_test, p_test, n_bins=8)
print("\nReliability plot data (a well-calibrated model has mp ≈ mt):")
print(f"  {'pred_prob':>10} {'true_frac':>10} {'count':>8}")
for a, b, c in zip(mp, mt, ct):
    print(f"  {a:>10.3f} {b:>10.3f} {c:>8d}")