Logistic Regression and Metrics

Access-change request REQ-10234 is sitting in the queue. May the workflow auto-process a permission change under policy line P-7, or must a human review it first?

The engine has to choose before the case is closed, so the true outcome isn't available yet. It needs a score that orders risk and a rule that turns that score into an action.

The previous chapter fitted latency from retrieved-context length with a design matrix, a loss, and a gradient-descent loop. We reuse that machinery, but now predict whether a request requires review. The score estimates a label; a separate policy chooses an action.

Logistic regression maps two evidence features to a score in (0, 1), then asks you to justify treating that score as a probability.^{[1]Reference 1The Elements of Statistical Learning.https://hastie.su.domains/ElemStatLearn/}

Eight constructed requests are enough to derive the gradient and check it in NumPy against scikit-learn.^{[2]Reference 2Array programming with NumPy.https://www.nature.com/articles/s41586-020-2649-2[3]Reference 3Scikit-learn: Machine Learning in Python.https://jmlr.org/papers/v12/pedregosa11a.html} A second slice, never used for fitting, is where a threshold's cost, the ranking, and the frequencies come apart. These rows teach the arithmetic. They aren't evidence for releasing a permission workflow.

From a latency number to a review decision

The previous target was continuous: "how many milliseconds will this evidence cost?" A residual of 2 ms still meant something. Now the label is binary: 1 means an access-change request needs human review before any permission change; 0 means policy evidence supports automatic handling.

For this fictional cost model, an incorrect automatic approval costs $120 in permission rollback and investigation. An unnecessary review costs $18 of queue time. Correct decisions have zero incremental cost. Those assumptions will choose the threshold later, not the loss used to fit the model. A real permission workflow must still enforce its authorization rules independently of this classifier.

Ordinary least squares can fit 0/1 labels as a linear probability model, but it has limitations:

1. Unbounded predictions: $\hat y=w^\top x+b$ can fall outside [0, 1]. Values such as -0.2 aren't valid probabilities.
2. Sensitivity to influential rows: unusual feature values and residuals can change the fitted boundary. The direction depends on the data; an extreme positive doesn't universally move a threshold rightward.
3. Conditional variance: for Bernoulli outcomes, $\operatorname{Var}(Y\mid x)=p(x)(1-p(x))$, which can vary with features. Constant variance isn't required to compute OLS or for unbiasedness under a correct conditional mean, but it is part of the classical BLUE guarantee and conventional standard-error formulas.

Squared error on bounded probabilities is also a legitimate scoring loss, the Brier score. Logistic regression gives us a Bernoulli likelihood with a linear logit; it isn't the only way to estimate probabilities.

A linear score can fall outside (0, 1), and a value such as -0.2 isn't a probability. Which smooth function keeps every score inside that interval and still has a derivative a training step can follow? The same split shows up in a guardrail or a router: the model emits a score, and a separate threshold turns that score into an action.

Squeeze the score into (0, 1)

Before choosing a threshold, we need a score that behaves like a probability: it must stay between 0 and 1, move upward when evidence favors review, and change smoothly enough for gradient descent. The sigmoid (or logistic) function does that:

Four properties matter for the later gradient:

• As $z \to +\infty$, $\sigma(z) \to 1$
• As $z \to -\infty$, $\sigma(z) \to 0$
• At $z = 0$, $\sigma(0) = 0.5$
• Its derivative is $\sigma'(z) = \sigma(z) \cdot (1 - \sigma(z))$, which is positive for finite $z$ and peaks at $0.25$ when $z = 0$. That tidy form is why the gradient of log-loss later collapses so cleanly.

Predict the shape before calculating. A large negative score should look almost like zero, z = 0 should land at one half, and a large positive score should approach one. Then compute five concrete values by hand:

import numpy as np

def sigmoid(z: np.ndarray) -> np.ndarray:
    return 1.0 / (1.0 + np.exp(-z))

scores = np.array([-4.0, -1.0, 0.0, 1.0, 4.0])
print("scores:", scores)
print("probabilities:", sigmoid(scores).round(3))

scores: [-4. -1.  0.  1.  4.]
probabilities: [0.018 0.269 0.5   0.731 0.982]

The code confirms the hand calculation. The curve is steep near zero and flat in both tails, which will matter when an update tries to move an already confident request.

View larger

The curve saturates quickly. A score of +4 produces p = 0.982; whether that number represents an honest frequency still needs validation on labels the fit did not see.

For finite real logits, sigmoid is strictly between 0 and 1 mathematically. Floating-point evaluation can round extreme results to 0 or 1. The first cell uses small logits for inspection; the trainer below uses a stable implementation and computes loss directly from logits.

The model remains linear inside the sigmoid:

The probability is $p = \sigma(z)$. The positive class is needs human review, so the sign and size of $z$ carry the evidence before the sigmoid bounds it.

The log-odds and logit link

The connection comes from odds, the ratio of the positive-label probability to the negative-label probability. An odds ratio would compare two odds, for example at two feature settings. Here the events are review required versus review not required, not the action the threshold will take:

When $p = 0.8$, the odds are $0.8 / 0.2 = 4$ (4 to 1 in favor of review). Taking the natural logarithm produces the logit, or log-odds:^{[1]Reference 1The Elements of Statistical Learning.https://hastie.su.domains/ElemStatLearn/}

Modeling the log-odds as a linear function $z = w^T x$ is the logit link function. Inverting that link gives the sigmoid activation directly:

At $z = 0$, the odds are 1:1 and $p = 0.5$. Adding one unit to $z$ multiplies the odds by $e \approx 2.718$.

The linear decision boundary hyperplane

Use $b=w_0$ and $\beta=[w_1,w_2]^\top$ for the feature coefficients, keeping the intercept separate. At $t=0.5$, review is selected when $z=b+\beta^\top x\ge0$. The equal-score boundary is

For $w_2\ne0$, we can write this line as:

For nonzero $\beta$, it is the normal vector pointing toward $p>0.5$. The signed orthogonal distance in these feature coordinates is:

Points sitting far from the hyperplane carry large $|z|$, driving $\sigma(z)$ into the saturated tails near 0 or 1. Changing the decision threshold to any $t \in (0, 1)$ shifts the decision boundary to the parallel hyperplane:

Lowering the threshold to $t = 0.20$ requires only $z \ge \log(0.25) \approx -1.386$, shifting the decision plane into negative logit space to route more borderline requests to human review.

If $w_2=0$ but $w_1\ne0$, the line is vertical. If $\beta=0$, the classifier is constant and has no separating hyperplane. Distance depends on feature units, so don't interpret it as a physical risk measure.

Eight historical requests

Each historical request has two scaled features produced from retrieved evidence. Read them as competing signals, not as finished decisions:

• x1 = ambiguity score: conflicting audit logs, missing requester context, or unclear scope description
• x2 = auto-policy support score: strength of evidence that policy line P-7 permits automatic handling

The label y = 1 when the completed case required human review before a permission change, and y = 0 when automatic handling was acceptable. The model must learn how those two signals combine.

These eight rows are constructed training data; the validation rows introduced later are separate:

The training fixture stays small so you can verify each number. Before seeing any weights, predict their signs: ambiguity should raise review risk, while policy support should lower it. R3 needed review despite only moderate ambiguity. R8 had substantial ambiguity, but strong policy support made automatic processing acceptable.

Add the intercept column of ones and build the design matrix X (8 × 3) and the label vector y.

Log-loss (binary cross-entropy) is the right loss for probabilities

If a request didn't need review (y = 0), a score of p = 0.93 gives the observed outcome only 1 - 0.93 = 0.07 probability. Its negative log is -log(0.07) ≈ 2.659. A score of p = 0.07 instead assigns that outcome probability 0.93, for a much smaller loss, -log(0.93) ≈ 0.073. A confident wrong prediction pays heavily.

Deriving binary cross-entropy from maximum likelihood

Each historical request has a binary label $y_i \in \{0, 1\}$. We model $y_i$ given features $x_i$ as a Bernoulli random trial with success probability $p_i = P(Y_i = 1 \mid x_i) = \sigma(w^T x_i)$. The Bernoulli probability mass function for a single example is:

When $y_i = 1$, $(1 - p_i)^0 = 1$ leaving $p_i$. When $y_i = 0$, $p_i^0 = 1$ leaving $1 - p_i$. Assuming conditional independence across the $N$ requests in our training set, the likelihood of observing the entire dataset is the joint product:

Multiplying hundreds of probabilities in $(0, 1)$ swiftly triggers floating-point underflow. Taking the natural logarithm converts the product into an additive sum of log-probabilities:

Maximizing this likelihood is equivalent to minimizing the Negative Log-Likelihood (NLL). Normalizing by $N$ gives the empirical Binary Cross-Entropy (BCE) loss:^{[4]Reference 4Pattern Recognition and Machine Learning.https://www.microsoft.com/en-us/research/publication/pattern-recognition-machine-learning/[5]Reference 5Machine Learning: A Probabilistic Perspective.https://probml.github.io/pml-book/book0.html}

For one hard label, define its target distribution $q_i=(1-y_i,y_i)$ and predicted distribution $r_i=(1-p_i,p_i)$. Then $H(q_i,r_i)=H(q_i)+D_{\mathrm{KL}}(q_i\parallel r_i)$ and $H(q_i)=0$ because this per-row target is a point mass. That doesn't make the population label distribution deterministic: our balanced training labels have nonzero marginal entropy. A sampled label isn't the full true conditional distribution. Expected log loss is minimized by that true conditional probability when the model can represent it.

Why not reuse Mean Squared Error?

We can fit with squared error, but composing it with sigmoid changes the optimization:

1. Nonconvexity in weights: squared error after sigmoid is generally nonconvex. BCE with a linear logit is convex. Its Hessian is $X^\top\operatorname{diag}(p_i(1-p_i))X/N$, positive definite at finite weights when $X$ has full column rank; redundant columns instead give flat directions.
2. Small gradients on confident mistakes: for squared error $(p-y)^2$, the logit derivative is $2(p-y)p(1-p)$. For half-squared error it is $(p-y)p(1-p)$. At $y=1,z=-8$, the former is about $-0.000670$, while BCE's derivative is about $-0.999665$. BCE preserves a strong learning signal on this confidently wrong prediction.

Strict convexity permits at most one finite minimizer. With complete or quasi-complete class separation, an unregularized optimum can require coefficients diverging to infinity. On perfectly separable data, loss can approach zero without attaining it at finite weights. Regularization can constrain that growth. Future accuracy still requires evaluation on representative unseen data.

Work one row by hand. For R4 (x1=3.5, x2=6, y=1) with a starting guess w = [0, 0, 0]:

• z = 0 + 0·3.5 + 0·6 = 0
• p = σ(0) = 0.5
• L = -[1 * log(0.5) + 0] = -log(0.5) ≈ 0.693

Now suppose the model improves in the correct direction and predicts $p = 0.82$ for the same row:

• L = -log(0.82) ≈ 0.198
• The loss dropped because the model became more confident in the correct direction, which is the behavior you want.

For code, compute the same loss from logits rather than taking log(sigmoid(z)). The expression log(1 + exp(z)) - y*z is the same loss, but evaluating exp(z) directly can overflow. NumPy's logaddexp(0, z) - y*z evaluates that softplus term stably even when z is very large in either direction.

The next check keeps the arithmetic in logit space. Predict which row should be cheap, which should be expensive, and whether extreme logits remain finite.

import numpy as np

def log_loss_from_logits(y: np.ndarray, z: np.ndarray) -> float:
    return float(np.mean(np.logaddexp(0.0, z) - y * z))

print("correct positive:", round(log_loss_from_logits(np.array([1.0]), np.array([2.0])), 3))
print("confident wrong:", round(log_loss_from_logits(np.array([1.0]), np.array([-4.0])), 3))
extreme = log_loss_from_logits(np.array([0.0, 1.0]), np.array([1000.0, -1000.0]))
print("extreme logits finite:", np.isfinite(extreme))

correct positive: 0.127
confident wrong: 4.018
extreme logits finite: True

The log-loss gradient collapses to (p - y) · x

The loss tells us how wrong a prediction is. The gradient tells each weight how to change the score. For R4, the starting probability error is 0.5 - 1 = -0.5. Multiplying by [1, 3.5, 6] gives [-0.5, -1.75, -3.0]. Subtracting this gradient raises the next score for that row. Why does this simple multiplication give the derivative?

Chain rule step by step

We have $L(p, y)$ with $p = \sigma(z)$ and linear score $z = w^T x = \sum_k w_k x_k$:

1. Loss with respect to probability: $$\frac{\partial L}{\partial p} = -\frac{y}{p} + \frac{1 - y}{1 - p} = \frac{-y(1 - p) + (1 - y)p}{p(1 - p)} = \frac{p - y}{p(1 - p)}$$
2. Probability with respect to score: $$\frac{\partial p}{\partial z} = \sigma'(z) = \sigma(z)(1 - \sigma(z)) = p(1 - p)$$
3. Score with respect to weight: $$\frac{\partial z}{\partial w_j} = x_j$$

Multiplying the first two terms produces a cancellation: the sigmoid derivative in the numerator eliminates the $p(1 - p)$ denominator from the cross-entropy derivative:

Compare the two label cases directly:

Both cases yield $\partial L / \partial z = p - y$. Combining this with $\partial z / \partial w_j = x_j$ gives the single-sample weight gradient:

Across the entire dataset of $N$ samples, stacking rows into design matrix $X \in \mathbb{R}^{N \times d}$, predicted probabilities $\hat{y} = \sigma(Xw) \in \mathbb{R}^N$, and label vector $y \in \mathbb{R}^N$ gives the batch gradient:^{[4]Reference 4Pattern Recognition and Machine Learning.https://www.microsoft.com/en-us/research/publication/pattern-recognition-machine-learning/[5]Reference 5Machine Learning: A Probabilistic Perspective.https://probml.github.io/pml-book/book0.html}

The bridge to linear regression

Notice the mathematical unity between this classifier and the linear regression model:

• Linear regression (MSE): $\nabla_w L = \frac{2}{N} X^T (\hat{y}_{\text{lin}} - y)$, where $\hat{y}_{\text{lin}} = Xw$
• Logistic regression (BCE): $\nabla_w L = \frac{1}{N} X^T (\hat{y}_{\text{log}} - y)$, where $\hat{y}_{\text{log}} = \sigma(Xw)$

Both use features transposed against prediction errors. The factor of 2 differs because the previous lesson used MSE; using half-MSE would make the prefactors match. Loss normalization changes gradient magnitudes and therefore the learning rate needed for the same update.

The logit is the Bernoulli canonical link in generalized linear models. The score-error term is simple, but weight updates still depend on feature values, averaging, step size, and any penalties or sample weights.

Plug in numbers for the starting point on R4 (p=0.5, y=1, x=[1, 3.5, 6]):

• $\text{grad}_{w0} = (0.5 - 1) \cdot 1 = -0.5$
• $\text{grad}_{w1} = (0.5 - 1) \cdot 3.5 = -1.75$
• $\text{grad}_{w2} = (0.5 - 1) \cdot 6 = -3.0$

The gradient points uphill, toward increasing loss. Gradient descent subtracts it, scaled by the learning rate, to move downhill.

The sign is the learning signal. If $y=1$ and $p=0.5$, $p-y$ is negative, so subtracting the gradient increases weights attached to positive features and raises the future score for similar requests. If $y=0$ and $p=0.8$, $p-y$ is positive, so subtracting the gradient lowers those weights and reduces the score.

Repeat that arithmetic for each row and average the gradients. You now have one logistic-regression update. The code keeps the linear-regression loop and swaps in the sigmoid prediction and (p - y) * x gradient.

import numpy as np

def sigmoid(z: float) -> float:
    return float(1.0 / (1.0 + np.exp(-z)))

def loss(w: np.ndarray, x: np.ndarray, y: float) -> float:
    z = x @ w
    return float(np.logaddexp(0.0, z) - y * z)

x = np.array([1.0, 3.5, 6.0])
y = 1.0
w = np.zeros(3)
analytic = (sigmoid(x @ w) - y) * x
eps = 1e-6
numeric = np.array([
    (loss(w + eps * np.eye(3)[j], x, y) - loss(w - eps * np.eye(3)[j], x, y)) / (2 * eps)
    for j in range(3)
])
print("analytic:", analytic.round(3))
print("numeric :", numeric.round(3))
print("match:", np.allclose(analytic, numeric, atol=1e-6))

analytic: [-0.5  -1.75 -3.  ]
numeric : [-0.5  -1.75 -3.  ]
match: True

View larger

Fit the eight rows in NumPy

Start with the full eight-row fit. The loop is the gradient-descent trainer from the latency chapter, with a sigmoid prediction and a (p - y) * x gradient. Before running it, predict the first update: the positive feature weights should rise from zero, while the intercept should stay near zero because the batch starts with equal positive and negative pressure.

import numpy as np

def sigmoid(z: np.ndarray) -> np.ndarray:
    # exp(-log(1 + exp(-z))) avoids overflow without clipping the score.
    return np.exp(-np.logaddexp(0.0, -z))

def log_loss_from_logits(y_true: np.ndarray, z: np.ndarray) -> float:
    if z.ndim != 1 or y_true.shape != z.shape:
        raise ValueError("labels and logits must both have shape (n,)")
    return float(np.mean(np.logaddexp(0.0, z) - y_true * z))

def fit_logistic(X: np.ndarray, y: np.ndarray, lr: float = 0.2, epochs: int = 3000, verbose: bool = True) -> np.ndarray:
    if X.ndim != 2 or len(X) == 0 or y.shape != (len(X),):
        raise ValueError("X must have shape (n, d), and y must have shape (n,), with n > 0")
    if not np.isfinite(X).all() or not np.isin(y, [0, 1]).all():
        raise ValueError("features must be finite and labels must be 0 or 1")
    n, d = X.shape
    w = np.zeros(d)
    for epoch in range(epochs):
        z = X @ w
        p = sigmoid(z)
        grad = (1 / n) * (X.T @ (p - y))
        w -= lr * grad
        if verbose and epoch in (0, 200, 500, 1000, 2999):
            loss = log_loss_from_logits(y, X @ w)
            print(f"after_epoch={epoch + 1:4d}  loss={loss:.4f}  w={np.round(w, 3)}")
    return w

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

# Historical access-change requests: [ambiguity, auto-policy support]
X_raw = np.array([
    [1.0, 1.0], [1.5, 2.0], [2.0, 3.0], [3.5, 6.0],
    [4.0, 7.0], [4.5, 5.0], [2.8, 4.0], [3.2, 8.0],
])
y = np.array([0, 0, 1, 1, 1, 1, 0, 0])

X = np.c_[np.ones((len(X_raw), 1)), X_raw]   # add intercept column

w = fit_logistic(X, y, lr=0.2, epochs=3000)
p = predict_proba(X, w)
print("\nFinal weights (w0, w1, w2):", w.round(3))
print("Probabilities:", p.round(3))
print("weights_match_reference", np.allclose(w, [-4.894, 3.399, -0.959], atol=0.01))

after_epoch=   1  loss=0.6830  w=[0.    0.069 0.075]
after_epoch= 201  loss=0.4648  w=[-2.205  1.615 -0.462]
after_epoch= 501  loss=0.4206  w=[-3.524  2.452 -0.686]
after_epoch=1001  loss=0.4092  w=[-4.356  3.023 -0.85 ]
after_epoch=3000  loss=0.4072  w=[-4.894  3.399 -0.959]

Final weights (w0, w1, w2): [-4.894  3.399 -0.959]
Probabilities: [0.079 0.153 0.274 0.777 0.88  0.996 0.687 0.156]
weights_match_reference True

Keep the label vector one-dimensional: X is (8, 3), w is (3,), and both X @ w and y are (8,). If an import produces y with shape (8, 1), multiplying it by (8,) logits broadcasts to (8, 8), pairing every label with every score. Averaging that matrix returns a scalar, but it isn't the rowwise log loss. The checks reject this mistake; convert a known single target column to y[:, 0] before fitting. Nonfinite features and labels outside {0, 1} also violate this trainer's Bernoulli model.

The negative intercept starts low when both signals are absent. Ambiguity's positive coefficient raises review risk, while the negative policy-support coefficient lowers it, holding the other feature fixed. Both feature weights rose on the first step because both features were larger among positive rows on average. Later steps separate their conditional effects; an initial gradient sign isn't the final coefficient sign.

Each progress line reports loss and weights after the same update, so snapshots stay comparable. The loss falls while the weights settle near the reference vector.

At $t=0.50$, R3 needs review but scores 0.274, while R7 doesn't but scores 0.687. A straight boundary can't create an isolated pocket. In real data, also inspect labels and missing features before interpreting errors as proof that the relationship needs a nonlinear model.

View larger

The fitted weights explain the training table. They don't pick a routing policy; that job moves to requests the optimizer never saw.

Select a threshold on held-out requests

The fitted model emits a score in (0, 1). A routing policy needs a threshold t: send the request to review when p >= t. Selecting t on the training rows would reward overfitting.

Use a validation slice with cases the optimizer never saw. Here a false negative means automatic processing when review was required and costs $120. A false positive means unnecessary human review and costs $18.

The code below tests three cuts on the same eight held-out rows. It counts caught reviews (tp), unnecessary reviews (fp), and missed reviews (fn). Precision measures the fraction of the review queue that truly needed review; recall measures the fraction of required reviews caught. Before running it, expect the lowest cut to catch more required reviews and create more queue work.

import numpy as np

def sigmoid(z: np.ndarray) -> np.ndarray:
    return 1.0 / (1.0 + np.exp(-z))

def metrics(y: np.ndarray, p: np.ndarray, threshold: float) -> tuple[float, float, float, int]:
    pred = (p >= threshold).astype(int)
    tp = int(((pred == 1) & (y == 1)).sum())
    fp = int(((pred == 1) & (y == 0)).sum())
    fn = int(((pred == 0) & (y == 1)).sum())
    precision = tp / (tp + fp) if tp + fp else 0.0
    recall = tp / (tp + fn) if tp + fn else 0.0
    f1 = 2 * precision * recall / (precision + recall) if precision + recall else 0.0
    cost = 120 * fn + 18 * fp
    return precision, recall, f1, cost

w = np.array([-4.894, 3.399, -0.959])
X_val_raw = np.array([
    [1.2, 2.0], [2.1, 4.0], [2.4, 3.0], [3.0, 5.0],
    [3.6, 4.0], [3.2, 7.0], [4.2, 5.0], [2.8, 6.0],
])
y_val = np.array([0, 0, 1, 0, 1, 0, 1, 1])
X_val = np.c_[np.ones(len(X_val_raw)), X_val_raw]
p_val = sigmoid(X_val @ w)

for t in (0.20, 0.50, 0.80):
    precision, recall, f1, cost = metrics(y_val, p_val, t)
    print(f"t={t:.2f} precision={precision:.3f} recall={recall:.3f} f1={f1:.3f} cost=${cost}")

print("p_val:", np.round(p_val, 3))
print("cheapest tested cut:", min((0.20, 0.50, 0.80), key=lambda t: metrics(y_val, p_val, t)[3]))

t=0.20 precision=0.667 recall=1.000 f1=0.800 cost=$36
t=0.50 precision=0.750 recall=0.750 f1=0.750 cost=$138
t=0.80 precision=1.000 recall=0.500 f1=0.667 cost=$240
p_val: [0.061 0.169 0.595 0.624 0.971 0.325 0.99  0.244]
cheapest tested cut: 0.2

Among these three cuts, t = 0.20 costs least. Compared with 0.50, it prevents one missed review and adds one unnecessary review: a net saving of $120 - $18 = $102. Once you select the threshold, freeze it and evaluate on an untouched test set. Validation cost was used to choose the policy, so it isn't an independent estimate of the selected policy's performance.^{[6]Reference 6Tuning the decision threshold for class predictionhttps://scikit-learn.org/stable/modules/classification_threshold.html}

There is also a probability-based decision rule, if the scores are calibrated. At p = 0.20, automatic processing has expected error cost 0.20 × $120 = $24; review has expected error cost 0.80 × $18 = $14.40. In general, review is cheaper when $18(1-p) \leq 120p$, or $p \geq 18/138 \approx 0.130$. That theoretical cut follows from our zero-cost-correct-decisions assumption. It needn't minimize observed cost on eight noisy cases, and it isn't justified by uncalibrated scores. For a live policy, estimate the full cost table and check reviewer capacity.

Precision, recall, F1, and class imbalance

For any threshold, fill the 2x2 confusion matrix. At t = 0.20, ask one concrete question for each validation request: did the route match the completed case?

View larger

At this cut the four cells represent concrete operational outcomes:

• True Positive (TP = 4): required reviews (V8, V3, V5, V7) routed to review. Correct labels don't prove that authorization or human review succeeds.
• False Positive (FP = 2): requests not requiring review (V6, V4) enter the queue, costing $18 each under our assumptions.
• True Negative (TN = 2): requests not requiring review (V1, V2) take the automatic route, with zero incremental cost in this toy table.
• False Negative (FN = 0): a request requiring review would take the automatic route, costing $120 in this toy model. This label alone doesn't establish that a breach occurred.

Six requests enter the queue, four of which need review: precision is 4/6 = 2/3. All four required reviews are caught: recall is 4/4 = 1. Their F1 is 2 × (2/3) × 1 / (2/3 + 1) = 0.8.

Read the matrix in the direction of the action:

• Precision = TP / (TP + FP): among requests routed to review, what fraction required it? Precision measures queue purity and alert trustworthiness.
• Recall = TP / (TP + FN): among requests requiring review, what fraction reached it? Recall measures safety coverage and catch rate.
• F1 = 2 * precision * recall / (precision + recall): the harmonic mean of precision and recall.

Precision is undefined when the selected queue is empty; recall is undefined when there are no positive labels. Our threshold helper uses zero for those summaries, a convention rather than an observed rate. Keep the class and queue counts alongside metrics. The reciprocal harmonic-mean expression below assumes positive precision and recall; use the count-based form or a declared zero-division convention at the boundaries.

Why F1 uses the harmonic mean

Why use the harmonic mean rather than the arithmetic average $(P + R) / 2$?

Consider a degenerate detector that flags only a single, glaringly obvious request where ambiguity is extreme ($x_1 = 10, x_2 = 0$). If that one request is positive, precision is $1 / 1 = 1.0$. But across 100 actual risky requests, recall is only $1 / 100 = 0.01$.

The arithmetic mean would be $(1.0 + 0.01) / 2 = 0.505$, making an utterly useless detector look moderately successful. The harmonic mean inverts this:

Because the harmonic mean is pulled strongly toward the smaller value, it plummets whenever either precision or recall fails. A high F1 requires both rates to remain high.

The accuracy trap under class imbalance

The four cells also explain why accuracy alone is dangerous when positive cases are rare. Try the extreme policy below: send every request to automatic processing.

import numpy as np

y_true = np.array([1] * 4 + [0] * 96)
predict_auto_process_everything = np.zeros(100, dtype=int)
accuracy = (predict_auto_process_everything == y_true).mean()
recall = ((predict_auto_process_everything == 1) & (y_true == 1)).sum() / (y_true == 1).sum()
print("accuracy:", round(float(accuracy), 3))
print("review recall:", round(float(recall), 3))

accuracy: 0.96
review recall: 0.0

The policy gets 96 of 100 rows right by exploiting the negative majority, yet it catches none of the four requests that required review. Accuracy answers a different question from review recall. When positive events are rare, a trivial majority-class baseline scores high accuracy while failing completely at its job.

ROC and precision-recall curves measure ranking

To draw an ROC (Receiver Operating Characteristic) curve, vary the threshold from high to low and plot recall (true-positive rate) on the y-axis against false-positive rate, FP / (FP + TN), on the x-axis. At t = 0.20, those coordinates are (2/4, 4/4) = (0.5, 1.0). Both classes must be present to estimate both rates.

AUC answers a ranking question: how often does a random positive score rank above a random negative score? Ties get half credit. AUC can tell you whether ordering is useful, but it never chooses a business threshold for you.

Before running the pairwise check, count the combinations: four positive and four negative validation rows should produce 16 comparisons.

import numpy as np

p_val = np.array([0.061, 0.169, 0.595, 0.624, 0.971, 0.325, 0.990, 0.244])
y_val = np.array([0, 0, 1, 0, 1, 0, 1, 1])
positive = p_val[y_val == 1]
negative = p_val[y_val == 0]
pair_scores = [(pos > neg) + 0.5 * (pos == neg) for pos in positive for neg in negative]
auc = float(np.mean(pair_scores))
print("positive-negative pairs:", len(pair_scores))
print("validation AUC:", round(auc, 3))

positive-negative pairs: 16
validation AUC: 0.812

This slice gives 13 ranking wins out of 16. The ordering is useful, even though one fixed threshold can still make an expensive mistake.

View larger

ROC-AUC versus PR-AUC under class imbalance

ROC curves measure threshold-independent discriminative ranking. The two axes normalize exclusively within each ground-truth class:

Population ROC rates depend on the class-conditional score distributions, not directly on the proportion of positives. Changing prevalence while preserving those distributions leaves the population ROC unchanged. Arbitrarily adding benign requests doesn't guarantee this: their score distribution may differ, and finite-sample estimates also vary. Uniformly replicating every negative row preserves an empirical ROC exactly; adding only easy low-score negatives doesn't.

Test that distinction on our eight scores. No model is retrained. Compare an unchanged negative-score mix with 100 extra negatives whose scores are all zero:

import numpy as np

p = np.array([0.061, 0.169, 0.595, 0.624, 0.971, 0.325, 0.990, 0.244])
y = np.array([0, 0, 1, 0, 1, 0, 1, 1])
positive = p[y == 1]
negative = p[y == 0]

def audit(labels: np.ndarray, scores: np.ndarray) -> tuple[float, float, float]:
    pos = scores[labels == 1]
    neg = scores[labels == 0]
    pairs = (pos[:, None] > neg[None, :]) + 0.5 * (pos[:, None] == neg[None, :])
    selected = scores >= 0.20
    precision = float(np.mean(labels[selected]))
    return float(np.mean(neg >= 0.20)), float(np.mean(pairs)), precision

same_mix_y = np.r_[np.ones(len(positive), dtype=int), np.zeros(25 * len(negative), dtype=int)]
same_mix_p = np.r_[positive, np.repeat(negative, 25)]
easy_y = np.r_[y, np.zeros(100, dtype=int)]
easy_p = np.r_[p, np.zeros(100)]

for name, labels, scores in (
    ("original", y, p),
    ("same-negative-mix", same_mix_y, same_mix_p),
    ("extra-easy-negatives", easy_y, easy_p),
):
    fpr, auc, precision = audit(labels, scores)
    print(name, "fpr", round(fpr, 4), "auc", round(auc, 4), "precision", round(precision, 4))

original fpr 0.5 auc 0.8125 precision 0.6667
same-negative-mix fpr 0.5 auc 0.8125 precision 0.0741
extra-easy-negatives fpr 0.0192 auc 0.9928 precision 0.6667

The unchanged mix lowers precision by adding false positives in the same proportion. The easy-negative addition instead inflates AUC while leaving the selected queue unchanged. State how an evaluation set was sampled before interpreting its ranking metrics.

That invariance becomes a blind spot when positives are rare. A classifier with recall 0.9 and false-positive rate 0.1 selects 90 true positives and 10 false positives from 100 positives and 100 negatives: precision 90 / 100 = 0.90. With only 100 positives among 10,000 requests (a 1:99 ratio), the exact same rates produce 90 true positives and 990 false positives: precision collapses to 90 / (90 + 990) ≈ 0.083.

Under the assumed unchanged rates, the ROC point stays (0.1, 0.9), but 91.7% of the queue consists of false positives. A precision-recall curve exposes that burden.^{[1]Reference 1The Elements of Statistical Learning.https://hastie.su.domains/ElemStatLearn/} For rare events, report AP together with counts, recall, cost, and reviewer capacity at the chosen threshold. Neither ROC-AUC nor an area under a PR curve alone measures operating effort.

View larger

Average precision (AP) summarizes the precision-recall curve with a stepwise weighted average of precision at each positive-label recall increase. People sometimes label AP as PR-AUC, but trapezoidal area under the same curve can produce a different number, so name the statistic you compute.^{[7]Reference 7average_precision_scorehttps://scikit-learn.org/stable/modules/generated/sklearn.metrics.average_precision_score.html}

Here the four positive rows enter at ranks 1, 2, 4, 6, with precisions 1, 1, 3/4, 4/6. Each adds one quarter of the total recall, so AP is their mean, about 0.854. This shortcut assumes distinct scores. Tied scores must enter together at one threshold; use average_precision_score for that general case.^{[7]Reference 7average_precision_scorehttps://scikit-learn.org/stable/modules/generated/sklearn.metrics.average_precision_score.html}

import numpy as np

p_val = np.array([0.061, 0.169, 0.595, 0.624, 0.971, 0.325, 0.990, 0.244])
y_val = np.array([0, 0, 1, 0, 1, 0, 1, 1])

# For distinct scores, AP is the mean precision at each positive row's rank.
order = np.argsort(-p_val)
y_sorted = y_val[order]
tp_cum = np.cumsum(y_sorted)
fp_cum = np.cumsum(1 - y_sorted)
precision_at = tp_cum / (tp_cum + fp_cum)
ap = float(precision_at[y_sorted == 1].mean())
print("validation average precision:", round(ap, 3))
print("validation ROC-AUC:   ", 0.812)

validation average precision: 0.854
validation ROC-AUC:    0.812

On this balanced slice AP and ROC-AUC are close. When positives become rare, report both. ROC-AUC can stay high while AP falls with positive prevalence and with the quality of ranking near the top of the queue.

Calibration: do scores represent frequencies?

A model can rank access-change requests well while producing over-confident or under-confident numbers. If requests scored near 0.80 require review only half the time, showing 80% to an operator misstates risk.

Logistic regression is often described as "already a probability model" because it optimizes log loss. That interpretation depends on the linear log-odds assumption and the data.

Shrinkage (scikit-learn's default C = 1.0, an L2 penalty) pulls coefficients toward zero, which can narrow the score range and move predictions toward the base rate. Even this unregularized eight-row fit isn't a frequency guarantee. Audit calibration on held-out labels.^{[1]Reference 1The Elements of Statistical Learning.https://hastie.su.domains/ElemStatLearn/}

A reliability diagram compares average score with observed positive frequency per bin. Our lowest bin has three requests with mean score 0.158, and one actually needs review (1/3 ≈ 0.333). Its absolute gap is about 0.175; because it contains 3/8 of the data, it contributes about 0.066 to expected calibration error (ECE).

We use binary positive-class ECE: bin the review probabilities, not the confidence in whichever class won. With $B_b$ denoting the rows in bin $b$, sum these weighted gaps:

ECE depends on sample size, prevalence, and bin choices. It has no universal pass threshold. Pair it with a Brier score, the mean squared difference between predicted probabilities and binary outcomes, when you need a bin-independent probability-quality summary. Brier score measures more than calibration: discrimination and outcome uncertainty matter too, so a lower Brier score doesn't by itself prove better calibration.^{[8]Reference 8Probability calibrationhttps://scikit-learn.org/stable/modules/calibration.html}

State the binning scheme and measure on held-out labels. Fit a post-hoc map on data separate from the classifier's fitting data, then verify the map on fresh labels. Platt scaling fits a one-dimensional logistic curve; isotonic regression fits a nondecreasing step function; temperature scaling divides logits by a learned positive scalar.^{[9]Reference 9On Calibration of Modern Neural Networkshttps://arxiv.org/abs/1706.04599} A strictly increasing map preserves binary ranking, while isotonic can introduce ties and change AUC.^{[8]Reference 8Probability calibrationhttps://scikit-learn.org/stable/modules/calibration.html} Reassess the numerical operating threshold on suitable validation data after calibration; don't tune it on the final test labels. Positive temperature scaling alone leaves the binary 0.5 decision boundary unchanged.

The next output makes the ECE bookkeeping visible rather than hiding it behind one number.

import numpy as np

p_val = np.array([0.061, 0.169, 0.595, 0.624, 0.971, 0.325, 0.990, 0.244])
y_val = np.array([0, 0, 1, 0, 1, 0, 1, 1])
edges = np.linspace(0.0, 1.0, 5)
ece = 0.0
for lo, hi in zip(edges[:-1], edges[1:]):
    mask = (p_val >= lo) & ((p_val < hi) if hi < 1.0 else (p_val <= hi))
    if mask.any():
        avg_p = p_val[mask].mean()
        frac_pos = y_val[mask].mean()
        ece += mask.mean() * abs(avg_p - frac_pos)
        closing = "]" if hi == 1.0 else ")"
        print(f"[{lo:.2f}, {hi:.2f}{closing} n={mask.sum()} avg_p={avg_p:.3f} frac_pos={frac_pos:.3f}")
print("validation ECE:", round(float(ece), 3))
print("validation Brier:", round(float(np.mean((p_val - y_val)**2)), 3))

[0.00, 0.25) n=3 avg_p=0.158 frac_pos=0.333
[0.25, 0.50) n=1 avg_p=0.325 frac_pos=0.000
[0.50, 0.75) n=2 avg_p=0.609 frac_pos=0.500
[0.75, 1.00] n=2 avg_p=0.980 frac_pos=1.000
validation ECE: 0.139
validation Brier: 0.158

In this sample, the first bin's observed review frequency exceeds its mean prediction (0.333 versus 0.158); the second is lower (0 versus 0.325). Those tiny bin counts don't establish population miscalibration. Calling the first bin "under-confident" would also obscure its high predicted confidence in the negative class. ECE uses an absolute gap regardless of direction.

Eight validation requests are enough for arithmetic, not a release gate. To see why, flip one disputed label and measure both metrics:

import numpy as np

def f1_at_half(y: np.ndarray, p: np.ndarray) -> float:
    pred = p >= 0.5
    tp = ((pred == 1) & (y == 1)).sum()
    fp = ((pred == 1) & (y == 0)).sum()
    fn = ((pred == 0) & (y == 1)).sum()
    return float(2 * tp / (2 * tp + fp + fn))

def ece(y: np.ndarray, p: np.ndarray) -> float:
    edges = np.linspace(0.0, 1.0, 5)
    total = 0.0
    for lo, hi in zip(edges[:-1], edges[1:]):
        mask = (p >= lo) & ((p < hi) if hi < 1.0 else (p <= hi))
        if mask.any():
            total += mask.mean() * abs(p[mask].mean() - y[mask].mean())
    return float(total)

p = np.array([0.061, 0.169, 0.595, 0.624, 0.971, 0.325, 0.990, 0.244])
base = np.array([0, 0, 1, 0, 1, 0, 1, 1])
one_flip = base.copy()
one_flip[3] = 1
for name, labels in (("base", base), ("one flip", one_flip)):
    print(f"{name:8s} f1={f1_at_half(labels, p):.3f} ece={ece(labels, p):.3f}")

base     f1=0.750 ece=0.139
one flip f1=0.889 ece=0.209

One label moves both results sharply. Treat this slice as a teaching audit, not a release gate. Two questions remain: what does each threshold cost, and do the scores describe frequencies? First, compare routing cost across thresholds:

View larger

Cost chooses an action. It doesn't say whether 0.87 means 87%. The reliability diagram asks that second question:

View larger

Verify the NumPy fit with scikit-learn

A from-scratch implementation gains a useful check by matching a maintained implementation under equivalent settings.^{[3]Reference 3Scikit-learn: Machine Learning in Python.https://jmlr.org/papers/v12/pedregosa11a.html} The September 21, 2026 library check uses NumPy 2.5.3 and scikit-learn 1.9.1. LogisticRegression defaults to L2 regularization (C=1.0, l1_ratio=0). Since 1.8, C=np.inf selects an unpenalized fit and penalty is deprecated; use C and l1_ratio with a compatible solver.^{[10]Reference 10LogisticRegressionhttps://scikit-learn.org/stable/modules/generated/sklearn.linear_model.LogisticRegression.html}

Set C=np.inf here so the library solve is comparable to our unregularized gradient-descent fit. The solver is L-BFGS rather than full-batch gradient descent, so the weights should be close, not bit-identical.

import numpy as np
from sklearn.linear_model import LogisticRegression
from sklearn.metrics import f1_score, roc_auc_score

X_raw = np.array([
    [1.0, 1.0], [1.5, 2.0], [2.0, 3.0], [3.5, 6.0],
    [4.0, 7.0], [4.5, 5.0], [2.8, 4.0], [3.2, 8.0],
])
y = np.array([0, 0, 1, 1, 1, 1, 0, 0])
X = np.c_[np.ones((len(X_raw), 1)), X_raw]
w_ours = np.array([-4.894, 3.399, -0.959])
X_val_raw = np.array([
    [1.2, 2.0], [2.1, 4.0], [2.4, 3.0], [3.0, 5.0],
    [3.6, 4.0], [3.2, 7.0], [4.2, 5.0], [2.8, 6.0],
])
y_val = np.array([0, 0, 1, 0, 1, 0, 1, 1])
X_val = np.c_[np.ones(len(X_val_raw)), X_val_raw]

clf = LogisticRegression(C=np.inf, solver="lbfgs", max_iter=5000, fit_intercept=False)
clf.fit(X, y)
w_sk = clf.coef_[0]
print("sklearn weights:", np.round(w_sk, 3))
print("our weights   :", w_ours)
print("weights close :", np.allclose(w_sk, w_ours, atol=0.05))
p_val = clf.predict_proba(X_val)[:, 1]
print("validation F1 @ 0.20:", round(f1_score(y_val, (p_val >= 0.20).astype(int)), 3))
print("validation AUC:", round(roc_auc_score(y_val, p_val), 3))

sklearn weights: [-4.914  3.414 -0.964]
our weights   : [-4.894  3.399 -0.959]
weights close : True
validation F1 @ 0.20: 0.8
validation AUC: 0.812

Close weights and matching metrics support our implementation on this fixture; they aren't a proof for other inputs. Compare objectives, logits, rank, and stopping tolerances when a match fails. The validation metrics describe eight constructed rows, not release performance.

Multiclass routing: the softmax extension

Binary logistic regression decides between two actions. Live policy engines often require more granular dispatch. Suppose our access-control engine sorts requests into three tiers:

• Class 0: Auto-Approve (standard permissions under established policy)
• Class 1: Standard Review (moderate risk ticket queue)
• Class 2: Security Escalation (severe risk requiring immediate review and credential freeze)

These are classifier labels. Permission checks and any credential actions are separately authorized workflow steps.

Linear heads and the softmax function

Instead of a single scalar logit, we assign each of the $K$ classes its own weight vector $w_k \in \mathbb{R}^d$ (including intercept):

In matrix form, the logit vector for one request is $z = W^T x$, where $W \in \mathbb{R}^{d \times K}$.

To convert unconstrained real logits $z \in \mathbb{R}^K$ into a valid probability distribution where each $p_k > 0$ and $\sum_{k=0}^{K-1} p_k = 1$, we apply the softmax function:^{[4]Reference 4Pattern Recognition and Machine Learning.https://www.microsoft.com/en-us/research/publication/pattern-recognition-machine-learning/[11]Reference 11Deep Learninghttps://www.deeplearningbook.org/}

Notice that when $K = 2$, softmax simplifies directly back to the binary sigmoid. Setting class logits to $z_0$ and $z_1$:

Binary logistic regression is two-class softmax regression with the reference class logit $z_0$ pinned to zero.

Categorical cross-entropy and its gradient

For a one-hot ground-truth vector $y \in \{0, 1\}^K$ ($y_k = 1$ for the true class, $0$ otherwise), the likelihood of a single observation under the categorical distribution is $\prod_{k=0}^{K-1} p_k^{y_k}$. Taking the negative log-likelihood across $N$ independent requests yields the Categorical Cross-Entropy loss:

Differentiating categorical cross-entropy with respect to class weight vector $w_k$ gives:

Stacking all $K$ weight vectors into weight matrix $W \in \mathbb{R}^{d \times K}$ produces the batch gradient:

where $X \in \mathbb{R}^{N \times d}$, $\hat{Y} \in \mathbb{R}^{N \times K}$ contains the softmax probabilities, and $Y \in \mathbb{R}^{N \times K}$ is the one-hot target matrix. The exact same error-residual form, features dotted with predicted minus actual, extends directly from scalar linear regression to multiclass classification.

Adding the same coefficient vector to every class head adds an identical feature-dependent value to every logit and leaves softmax unchanged. Thus unrestricted class weights aren't identifiable from probabilities alone. A reference class or an appropriate constraint removes that redundancy; convexity doesn't imply a unique unrestricted weight matrix.

Numerical stability: the log-sum-exp trick

Direct exponentials overflow around logits 709.8 in float64 or 88.7 in float32. Subtracting $c=\max_jz_j$ preserves softmax mathematically and stabilizes those exponentials:

With representable finite differences, $z_k-c\le0$, so the largest exponential is 1. Tiny probabilities can still underflow or round to zero; invalid inputs and extreme subtraction overflow need separate handling. Compute cross-entropy from shifted logits, not clipped probabilities. Clipping would cap the penalty for an extremely wrong prediction.

Fit a three-class routing classifier on six historical requests using NumPy:

import numpy as np

def softmax(Z: np.ndarray) -> np.ndarray:
    shift_Z = Z - np.max(Z, axis=1, keepdims=True)
    exp_Z = np.exp(shift_Z)
    return exp_Z / np.sum(exp_Z, axis=1, keepdims=True)

def categorical_cross_entropy_from_logits(Y_onehot: np.ndarray, Z: np.ndarray) -> float:
    shifted = Z - np.max(Z, axis=1, keepdims=True)
    log_normalizer = np.log(np.sum(np.exp(shifted), axis=1, keepdims=True))
    log_probabilities = shifted - log_normalizer
    return float(-np.mean(np.sum(Y_onehot * log_probabilities, axis=1)))

# 6 requests: [ambiguity, auto-policy support] + intercept column
# Target routing classes: 0: Auto-Approve, 1: Standard Review, 2: Escalation
X_raw = np.array([
    [1.0, 8.0],
    [1.5, 7.0],
    [2.5, 4.0],
    [3.2, 5.0],
    [4.8, 1.5],
    [5.0, 1.0],
])
X = np.c_[np.ones(len(X_raw)), X_raw]
y_classes = np.array([0, 0, 1, 1, 2, 2])
K = 3
N, d = X.shape
Y_onehot = np.eye(K)[y_classes]

W = np.zeros((d, K))
lr = 0.2
for epoch in range(1000):
    P = softmax(X @ W)
    grad = (1 / N) * (X.T @ (P - Y_onehot))
    W -= lr * grad

P_final = softmax(X @ W)
preds = np.argmax(P_final, axis=1)
loss = categorical_cross_entropy_from_logits(Y_onehot, X @ W)
print("final loss:", round(loss, 4))
print("predicted classes:", preds.tolist())
print("all correct:", bool((preds == y_classes).all()))
print("sample probabilities row 0 (auto):", P_final[0].round(3).tolist())
print("sample probabilities row 4 (escalate):", P_final[4].round(3).tolist())

extreme_logits = np.array([[1000.0, 0.0, -1000.0]])
extreme_target = np.array([[0.0, 0.0, 1.0]])
extreme_loss = categorical_cross_entropy_from_logits(extreme_target, extreme_logits)
assert np.isfinite(extreme_loss) and np.isclose(extreme_loss, 2000.0)
print("extreme wrong log loss:", round(extreme_loss, 1))

final loss: 0.0055
predicted classes: [0, 0, 1, 1, 2, 2]
all correct: True
sample probabilities row 0 (auto): [1.0, 0.0, 0.0]
sample probabilities row 4 (escalate): [0.0, 0.005, 0.995]
extreme wrong log loss: 2000.0

"All correct" scores the six training rows. The printed [1.0, 0.0, 0.0] is rounded, not calibrated certainty. The extreme wrong prediction costs 2,000 nats; clipping its target probability at 1e-15 would misleadingly report only about 34.54. This separable toy fit also doesn't establish finite maximum-likelihood coefficients.

The logistic and softmax models give linear log-odds scores followed by normalized exponential curves. Before moving to trees, compare that objective with two other uses of feature geometry: a separating margin and a covariance kernel.

Margins, kernels, and predictive uncertainty

A support vector machine (SVM) can classify the same access requests with a linear score $f(x)=w^\top x+b$. Its question differs from logistic regression's: can the classes be separated with a wide geometric buffer around the boundary?

In the usual soft-margin SVM, support vectors lie on or inside the margin, or on the wrong side of the decision boundary. Rows comfortably beyond the margin don't contribute hinge-loss pressure.^{[12]Reference 12Support Vector Machineshttps://scikit-learn.org/stable/modules/svm.html}

Encode needs review as $y=+1$ and review not required as $y=-1$. The functional margin is $yf(x)$: positive means the request is classified correctly. Dividing by the feature-weight norm gives the signed geometric margin for a linear model.

A standard hinge loss charges:

A positive request with score 0.4 is classified correctly, but its margin is only 0.4, so hinge loss is 0.6. A positive score of 1.4 clears the unit margin and incurs zero hinge loss.

The logistic loss on the second request is about 0.220; hinge loss is zero after the unit functional margin. A usual soft-margin objective balances violations with weight regularization, for example $\frac12\|w\|^2+C\sum_i\ell_{\mathrm{hinge}}(y_i,f(x_i))$. LinearSVC defaults to squared hinge and also regularizes the intercept through liblinear. Name the loss, penalty, and normalization when comparing implementations.^{[12]Reference 12Support Vector Machineshttps://scikit-learn.org/stable/modules/svm.html}

A linear SVM still can't carve out an isolated exception that requires a curved boundary. A kernel method replaces an explicit feature expansion with a pairwise similarity function. The radial basis function (RBF) kernel is one common choice:

Here $\ell$ is the length scale. Identical points have similarity 1. With $\ell=1$, points one unit apart have similarity 0.607, and points three units apart have similarity 0.011.

An RBF SVM can draw nonlinear boundaries around feature neighborhoods without materializing an infinite-dimensional feature map. Its decision score still isn't a calibrated class probability. Apply a separate held-out calibration procedure before displaying percentages, and scale features before treating Euclidean distance as meaningful.^{[12]Reference 12Support Vector Machineshttps://scikit-learn.org/stable/modules/svm.html}

A Gaussian process (GP) can use that same RBF kernel for a different job: define covariance between unknown function values and return a predictive mean together with uncertainty. Switch briefly to a regression example, not the binary review label: observe a continuous value $y=2$ at input 0. Assume zero prior mean, unit prior variance, length scale 1, and Gaussian observation-noise variance 0.04. At the observed input, the posterior mean is 2/1.04 ≈ 1.923 and latent variance is 1 - 1/1.04 ≈ 0.0385. At a new input the same conditioning calculation gives:^{[13]Reference 13Gaussian Processes for Machine Learning: Chapter 2, Regressionhttps://gaussianprocess.org/gpml/chapters/RW2.pdf}

Near the observation, the mean stays close to 2 and uncertainty shrinks. Far away, similarity vanishes, the mean returns toward its zero prior, and standard deviation rises toward 1. This variance describes the latent function; predicting a new noisy observation adds 0.04 to it. GP classification requires a non-Gaussian likelihood and doesn't use these regression equations unchanged.

This is predictive uncertainty under the stated kernel, noise, and prior assumptions. Check calibration on shifted traffic before relying on the intervals. Exact dense GP regression also requires roughly $O(n^3)$ fitting time and $O(n^2)$ storage for $n$ observations, so it fits smaller data problems better than large unapproximated serving datasets.^{[13]Reference 13Gaussian Processes for Machine Learning: Chapter 2, Regressionhttps://gaussianprocess.org/gpml/chapters/RW2.pdf}

The following dependency-free comparison makes both distinctions concrete: it prints the two loss objectives and the one-observation GP posterior:

from math import exp, log1p, sqrt

examples = [
    ("required-review", 1, 0.4),
    ("confident-review", 1, 1.4),
    ("missed-review", 1, -0.2),
    ("safe-auto", -1, -1.2),
]

for name, label, score in examples:
    margin = label * score
    hinge = max(0.0, 1.0 - margin)
    logistic = log1p(exp(-margin))
    print(f"{name:16} margin={margin:4.1f} hinge={hinge:.3f} logistic={logistic:.3f}")

print("distance  kernel  gp_mean  gp_std")
for distance in [0.0, 1.0, 3.0]:
    similarity = exp(-(distance**2) / 2.0)
    mean = 2.0 * similarity / 1.04
    standard_deviation = sqrt(1.0 - similarity**2 / 1.04)
    print(f"{distance:8.1f}  {similarity:.3f}   {mean:.3f}    {standard_deviation:.3f}")

required-review  margin= 0.4 hinge=0.600 logistic=0.513
confident-review margin= 1.4 hinge=0.000 logistic=0.220
missed-review    margin=-0.2 hinge=1.200 logistic=0.798
safe-auto        margin= 1.2 hinge=0.000 logistic=0.263
distance  kernel  gp_mean  gp_std
     0.0  1.000   1.923    0.196
     1.0  0.607   1.166    0.804
     3.0  0.011   0.021    1.000

The margin rows compare two loss terms; fitted models also depend on their penalties. In our one-observation stationary-RBF example, uncertainty rises with distance from that observation. With other kernels, priors, and several observations, uncertainty needn't be a simple monotone function of distance.

If a run diverges from its expected numbers, diagnose the boundary between data, objective, and implementation before blaming the optimizer.

When the numbers look wrong

Try it yourself

Use the same training and validation rows. Follow the order: verify the update, stress the optimizer, choose a policy, then test how fragile the probability story is. Fit only on training; choose policy only on validation.

1. Start with w = [0, 0, 0] and compute the gradient vector on the full batch by hand for the initial step (lr = 0.2). Verify it matches the first printed line of the NumPy run: w = [0, 0.069, 0.075].
2. Change the learning rate to 2.0. Print the first ten losses instead of only the final loss. Does each update improve the objective? Explain any increases in terms of step size and feature scale.
3. Sweep thresholds [0.2, 0.35, 0.5, 0.65, 0.8] on validation and print precision, recall, F1, and cost. Which threshold minimizes the stated cost?
4. Bin validation probabilities into four equal-width bins, compute ECE by hand, and compare with code output. Repeat with eight bins and explain why the number changes.
5. Flip one validation label and measure F1 and ECE again. Write the argument for gathering a larger, later-in-time validation set before publishing probability claims.

Check your calculations: the initial mean gradient is [0, -0.34375, -0.375]. The five thresholds give costs $36, $138, $138, $240, $240, respectively. Four-bin ECE is about 0.139; eight-bin ECE is about 0.154 using the displayed three-decimal scores. The changed bins expose gaps that previously canceled within a bin. Neither value settles calibration on such a small sample.

Turn two evidence scores into a probability, choose a threshold from the stated dollar costs, and keep a ranking result separate from a calibration gap. On the training table, R3 still scores only 0.274 at the fitted weights: moderate ambiguity, moderate policy support, but it needed review. A linear score has no way to carve out that exception without also moving nearby automatic requests. Axis-aligned splits on this same table are how you check whether a rule can isolate R3.