Metrics — Regression & Classification

# Every regression metric from scratch on a toy fit (EQ V3.1-V3.4).
import numpy as np
rng = np.random.default_rng(0)

n = 60
x = np.linspace(0, 10, n)
y = 3.0 * x + 5.0 + rng.normal(0, 2.0, n)     # truth: a noisy line
# Fit y = a*x + b by least squares (this is what MSE training would find).
A = np.vstack([x, np.ones_like(x)]).T
a, b = np.linalg.lstsq(A, y, rcond=None)[0]
yhat = a * x + b

e = y - yhat                                   # residuals: the raw material
mse  = np.mean(e**2)
rmse = np.sqrt(mse)
mae  = np.mean(np.abs(e))
mape = 100 * np.mean(np.abs(e / y))            # y is safely far from 0 here
r2   = 1 - np.sum(e**2) / np.sum((y - y.mean())**2)

print(f"fitted line : y = {a:.2f}*x + {b:.2f}")
print(f"MSE  = {mse:.3f}")
print(f"RMSE = {rmse:.3f}   (same units as y)")
print(f"MAE  = {mae:.3f}    (<= RMSE, always)")
print(f"MAPE = {mape:.2f} %")
print(f"R2   = {r2:.3f}     (1 = perfect, 0 = no better than the mean)")
plot_scatter(x, y)                              # the data the line was fit to

# One outlier: MAE shrugs, RMSE leaps. And MAPE near zero detonates.
import numpy as np

e = np.ones(10)                                # ten residuals of size 1
print("clean residuals:", e.astype(int))
print(f"  MAE  = {np.mean(np.abs(e)):.3f}   RMSE = {np.sqrt(np.mean(e**2)):.3f}")

e_out = e.copy(); e_out[0] = 10                # turn ONE into a size-10 miss
print("\nwith one size-10 outlier:")
print(f"  MAE  = {np.mean(np.abs(e_out)):.3f}   RMSE = {np.sqrt(np.mean(e_out**2)):.3f}")
print(f"  RMSE rose {np.sqrt(np.mean(e_out**2))/np.sqrt(np.mean(e**2)):.2f}x; "
      f"MAE rose only {np.mean(np.abs(e_out))/np.mean(np.abs(e)):.2f}x")

# The MAPE trap: identical absolute errors, but one true value sits near zero.
y     = np.array([100., 100., 100.,   0.01])
yhat  = np.array([101., 101., 101.,   0.50])
print("\nMAPE per row (%):", np.round(100*np.abs((y-yhat)/y), 1))
print(f"overall MAPE = {100*np.mean(np.abs((y-yhat)/y)):.0f} %  "
      "<- one near-zero label hijacks the whole average")

# Confusion matrix -> precision, recall, F1, accuracy, all from scratch.
import numpy as np
rng = np.random.default_rng(2)

n = 1000
y = rng.integers(0, 2, n)                          # true labels
# scores: positives score higher on average, but the classes overlap
score = rng.normal(0, 1, n) + 1.3 * y
thr = 0.5
pred = (score > thr).astype(int)

TP = int(np.sum((pred == 1) & (y == 1)))
FP = int(np.sum((pred == 1) & (y == 0)))
FN = int(np.sum((pred == 0) & (y == 1)))
TN = int(np.sum((pred == 0) & (y == 0)))
print(f"confusion: TP={TP}  FP={FP}  FN={FN}  TN={TN}")

precision = TP / (TP + FP)                         # of predicted positives, how many real
recall    = TP / (TP + FN)                         # of real positives, how many caught
f1        = 2 * precision * recall / (precision + recall)
accuracy  = (TP + TN) / n
print(f"precision = {precision:.3f}")
print(f"recall    = {recall:.3f}")
print(f"F1        = {f1:.3f}   (harmonic mean of the two)")
print(f"accuracy  = {accuracy:.3f}")

# The accuracy paradox: 99% accurate and completely useless.
import numpy as np
rng = np.random.default_rng(5)

n = 10000
y = (rng.random(n) < 0.01).astype(int)            # 1% positive (rare disease)
print(f"true positives in data: {y.sum()} / {n}  ({100*y.mean():.1f}%)\n")

# "Model" A: predict the majority class (negative) for everyone.
predA = np.zeros(n, dtype=int)
accA = (predA == y).mean()
TP_A = int(np.sum((predA==1)&(y==1)))
print(f"ALWAYS-NEGATIVE:  accuracy = {accA:.3f}  but caught {TP_A} of {y.sum()} positives")
print("  precision/recall/F1 are 0/0 -> the metric that 'works' is a mirage.\n")

# "Model" B: a real but imperfect detector.
score = rng.normal(0, 1, n) + 2.5 * y
predB = (score > 1.5).astype(int)
TP=int(np.sum((predB==1)&(y==1))); FP=int(np.sum((predB==1)&(y==0)))
FN=int(np.sum((predB==0)&(y==1))); TN=int(np.sum((predB==0)&(y==0)))
prec = TP/(TP+FP) if TP+FP else 0
rec  = TP/(TP+FN) if TP+FN else 0
f1   = 2*prec*rec/(prec+rec) if prec+rec else 0
print(f"REAL DETECTOR:    accuracy = {(TP+TN)/n:.3f}")
print(f"  precision = {prec:.3f}  recall = {rec:.3f}  F1 = {f1:.3f}")
print("Same-ish accuracy, but only F1/precision/recall reveal which model works.")

# Log loss vs Brier: confidence is rewarded only when it's right (EQ V3.9).
import numpy as np

def log_loss(y, p):
    p = np.clip(p, 1e-12, 1 - 1e-12)              # guard ln(0) = -inf
    return -np.mean(y*np.log(p) + (1-y)*np.log(1-p))

def brier(y, p):
    return np.mean((p - y)**2)

y = np.array([1, 1, 0, 0])                        # two positives, two negatives
confident_right = np.array([0.95, 0.90, 0.05, 0.10])
hedging         = np.array([0.55, 0.55, 0.45, 0.45])
confident_wrong = np.array([0.05, 0.10, 0.95, 0.90])

for name, p in [("confident & right", confident_right),
                ("hedging (~0.5)   ", hedging),
                ("confident & WRONG", confident_wrong)]:
    print(f"{name}:  log loss = {log_loss(y,p):.3f}   Brier = {brier(y,p):.3f}")

print("\nlog loss of a single confident-correct 0.9  :", round(-np.log(0.9), 3))
print("log loss of a single coin-flip      0.5     :", round(-np.log(0.5), 3))
print("log loss explodes for confident errors; Brier stays bounded by 1.")