Metrics

This module scores a model after it has been fitted: the confusion counts and the rates built from them for binary classification, the area under the ROC curve, probabilistic scores, regression errors, and the silhouette for a clustering with no labels.

Conventions. Every score reads its arguments as (truth, prediction): y first, then p. For classification both are vectors of 0 and 1 of the same length, with 1 meaning positive. For probabilistic scores p holds predicted probabilities of the positive class. For regression y and p are real vectors. The loss functions in Prelude.Learning take (prediction, truth), the order an optimiser differentiates in; the scores here turn that around once so every call site reads the same way.

Division by zero is guarded. Where a denominator can be zero, such as precision for a classifier that predicts no positives, the zero is replaced by 1, so the score comes back as 0 rather than NaN. The classification precision and recall are named precisionScore and recallScore, because precision already means the inverse covariance matrix in Prelude.Multivariate.

open Prelude.Metrics

The four counts

A binary confusion matrix is four sums over the observations, one for each combination of truth and prediction:

\[\mathrm{TP} = \sum_i y_i\, p_i, \quad \mathrm{FP} = \sum_i (1 - y_i)\, p_i, \quad \mathrm{FN} = \sum_i y_i (1 - p_i), \quad \mathrm{TN} = \sum_i (1 - y_i)(1 - p_i)\]

hits

hits(y, p) · y true labels, 0/1 vector n · p predicted labels, 0/1 vector n → true positives scalar

Let \(y \in \{0, 1\}^n\) be the true labels and \(p \in \{0, 1\}^n\) the predicted labels. Then

\[\operatorname{hits}(y, p) = \mathrm{TP} = \sum_{i=1}^{n} y_i\, p_i .\]

The number of positives predicted as positive.

Intuition The product \(y_i p_i\) is 1 only when the truth and the prediction are both positive, so summing it counts the positives the model caught. It is the shared numerator of precision and recall.

Use it when you build a confusion matrix or a custom score from counts. Both arguments must be hard 0/1 labels: with probabilities in p the sum becomes a soft, expected count, so apply a threshold first. For rates rather than counts use precisionScore and recallScore.

classificationmodel-evaluation

def hits(y, p)       ≔ einsum("i,i->", y, p)

falseCalls

falseCalls(y, p) · y true labels vector n · p predicted labels vector n → false positives scalar

Let \(y \in \{0, 1\}^n\) be the true labels and \(p \in \{0, 1\}^n\) the predicted labels. Then

\[\operatorname{falseCalls}(y, p) = \mathrm{FP} = \sum_{i=1}^{n} (1 - y_i)\, p_i .\]

The number of negatives predicted as positive.

Intuition False alarms: the model said positive and the truth was negative. The factor \((1 - y_i)\) flips the truth, so the product is 1 only on a negative that was called positive.

Use it when false positives carry the cost you care about, such as a spam filter binning real mail, a fraud rule blocking good customers, or an alert waking an engineer. It feeds precisionScore, and the false positive rate is 1 − specificity(y, p). Watch it together with missed when you move a decision threshold, because lowering the threshold trades one for the other.

classificationerror-analysisthreshold-selection

def falseCalls(y, p) ≔ einsum("i,i->", 1 − y, p)

missed

missed(y, p) · y true labels vector n · p predicted labels vector n → false negatives scalar

Let \(y \in \{0, 1\}^n\) be the true labels and \(p \in \{0, 1\}^n\) the predicted labels. Then

\[\operatorname{missed}(y, p) = \mathrm{FN} = \sum_{i=1}^{n} y_i\, (1 - p_i) .\]

The number of positives predicted as negative.

Intuition The positives the model let through: the truth is 1 and the prediction is 0. The factor \((1 - p_i)\) flips the prediction, so only the misses survive the product.

Use it when a miss is the expensive error, such as a disease not detected, a defect shipped, or fraud not flagged. It is the count that recallScore rewards you for driving down. When tuning a threshold, lowering it usually turns misses into falseCalls, so track both.

classificationerror-analysisthreshold-selection

def missed(y, p)     ≔ einsum("i,i->", y, 1 − p)

rejected

rejected(y, p) · y true labels vector n · p predicted labels vector n → true negatives scalar

Let \(y \in \{0, 1\}^n\) be the true labels and \(p \in \{0, 1\}^n\) the predicted labels. Then

\[\operatorname{rejected}(y, p) = \mathrm{TN} = \sum_{i=1}^{n} (1 - y_i)(1 - p_i) .\]

The number of negatives predicted as negative.

Intuition The negatives correctly left alone: both flags are flipped, so the product is 1 only when truth and prediction are both 0. On imbalanced data this is usually by far the largest of the four counts.

Use it when you compute specificity, accuracy or matthews, or report the full confusion matrix. Because it is so large on rare-positive problems, a score that includes it (such as accuracy) can look excellent for a model that finds no positives at all. Check recallScore or matthews alongside it.

classificationmodel-evaluationclass-imbalance

def rejected(y, p)   ≔ einsum("i,i->", 1 − y, 1 − p)

Building blocks

called(y, p) — Let \(y, p \in \{0, 1\}^n\). Then \(\operatorname{called}(y, p) = \mathrm{TP} + \mathrm{FP} = \sum_i p_i\). The number of predicted positives.
actual(y) — Let \(y \in \{0, 1\}^n\). Then \(\operatorname{actual}(y) = \sum_{i=1}^{n} y_i\). The number of true positives in the data.

The scores

The classification rates are ratios of the four counts:

\[\text{accuracy} = \frac{\mathrm{TP} + \mathrm{TN}}{n}, \quad \text{precision} = \frac{\mathrm{TP}}{\mathrm{TP} + \mathrm{FP}}, \quad \text{recall} = \frac{\mathrm{TP}}{\mathrm{TP} + \mathrm{FN}}, \quad \text{specificity} = \frac{\mathrm{TN}}{\mathrm{TN} + \mathrm{FP}}\] \[F_\beta = \frac{(1 + \beta^2)\,P\,R}{\beta^2 P + R}, \qquad \mathrm{MCC} = \frac{\mathrm{TP}\cdot\mathrm{TN} - \mathrm{FP}\cdot\mathrm{FN}}{\sqrt{(\mathrm{TP}+\mathrm{FP})(\mathrm{TP}+\mathrm{FN})(\mathrm{TN}+\mathrm{FP})(\mathrm{TN}+\mathrm{FN})}}\]

\(F_1\) is the harmonic mean of precision and recall, so it cannot be high unless both are. The Matthews correlation coefficient uses all four counts, none of which can be ignored, so it stays informative on an imbalanced data set where accuracy does not.

accuracy

accuracy(y, p) · y true labels vector n · p predicted labels vector n → accuracy scalar

\[\operatorname{accuracy}(y, p) = \frac{\mathrm{TP} + \mathrm{TN}}{n} .\]

The share of observations predicted correctly.

Intuition The fraction of predictions that are right, counting both classes. It treats a false alarm and a miss as equally bad and is dominated by the majority class: always predicting “negative” on data with 1% positives scores 0.99.

Use it when the classes are roughly balanced and both errors cost about the same, as in a quick benchmark. On imbalanced data use matthews, f1 or auc. When one kind of error matters more, look at precisionScore and recallScore separately.

classificationmodel-evaluation

def accuracy(y, p)       ≔ (hits(y, p) + rejected(y, p)) / count(y)

precisionScore

precisionScore(y, p) · y true labels vector n · p predicted labels vector n → precision scalar

Let \(y \in \{0, 1\}^n\) be the true labels and \(p \in \{0, 1\}^n\) the predicted labels. Let \(\mathrm{TP}, \mathrm{FP}, \mathrm{FN}, \mathrm{TN}\) be the four counts of hits, falseCalls, missed and rejected, and let \(c = \mathrm{TP} + \mathrm{FP}\). Then

\[\operatorname{precisionScore}(y, p) = P = \frac{\mathrm{TP}}{c + [\, c = 0 \,]} .\]

The share of predicted positives that are truly positive. Returns 0 when nothing is predicted positive.

Intuition Of everything the model flagged, how much was really positive? It measures how far a positive call can be trusted: precision 0.9 means one flag in ten is a false alarm. It says nothing about the positives the model never flagged.

Use it when acting on a flag is expensive, as with manual review queues, blocked transactions or recommended items. A model can raise precision by flagging only its surest cases, so report it with recallScore, or combine the two with f1 or fbeta. The divisor is guarded, so a model that flags nothing scores 0 rather than NaN.

classificationmodel-evaluationthreshold-selection

def precisionScore(y, p) ≔ hits(y, p) / safe(called(y, p))

recallScore

recallScore(y, p) · y true labels vector n · p predicted labels vector n → recall scalar

\[\operatorname{recallScore}(y, p) = R = \frac{\mathrm{TP}}{a + [\, a = 0 \,]} .\]

The share of true positives that are found (sensitivity). Returns 0 when the data contain no positives.

Intuition Of all the real positives, how many did the model find? Also called sensitivity or the true positive rate. It ignores false alarms entirely, so flagging everything gives a recall of 1.

Use it when missing a positive is the costly error: screening tests, safety defects, or search where the user must see every relevant document. Pair it with precisionScore, or use fbeta with b above 1 to weight it more heavily. It is the vertical axis of an ROC curve. The divisor is guarded, so data with no positives give 0.

classificationmodel-evaluationthreshold-selection

def recallScore(y, p)    ≔ hits(y, p) / safe(actual(y))

specificity

specificity(y, p) · y true labels vector n · p predicted labels vector n → specificity scalar

\[\operatorname{specificity}(y, p) = \frac{\mathrm{TN}}{(n - a) + [\, n - a = 0 \,]} .\]

The share of true negatives that are predicted negative. Returns 0 when the data contain no negatives.

Intuition Of all the real negatives, how many did the model correctly leave alone? This is the true negative rate. One minus it is the false positive rate, the horizontal axis of an ROC curve.

Use it when false alarms among negatives matter, as in medical tests (reported next to sensitivity, which is recallScore) or quality checks that should not reject good parts. On heavily imbalanced data a high specificity is easy to get, so check precisionScore to see how many flags are real. The divisor is guarded, so data with no negatives give 0.

classificationmodel-evaluation

def specificity(y, p)    ≔ rejected(y, p) / safe(count(y) − actual(y))

f1

f1(y, p) · y true labels vector n · p predicted labels vector n → F₁ scalar

Let \(y \in \{0, 1\}^n\) be the true labels and \(p \in \{0, 1\}^n\) the predicted labels. Let \(P = \operatorname{precisionScore}(y, p)\) and \(R = \operatorname{recallScore}(y, p)\). Then

\[\operatorname{f1}(y, p) = \frac{2\, P\, R}{(P + R) + [\, P + R = 0 \,]} .\]

The harmonic mean of precision and recall. Returns 0 when both are 0.

Intuition The harmonic mean punishes imbalance between the two rates: precision 1.0 with recall 0.1 gives an \(F_1\) of about 0.18, where the arithmetic mean would say 0.55. It is high only when both are high, and it never looks at true negatives.

Use it when positives are rare and you want one number balancing false alarms and misses, as in information extraction, entity matching or rare-event detection. If one error matters more than the other, use fbeta. Because true negatives are ignored, the score depends on which class you call positive; matthews does not. It is guarded to 0 when precision and recall are both 0.

classificationmodel-evaluationclass-imbalance

def f1(y, p)       ≔ 2 · precisionScore(y, p) · recallScore(y, p) / safe(precisionScore(y, p) + recallScore(y, p))

fbeta

fbeta(y, p, b) · y true labels vector n · p predicted labels vector n · b β, the weight of recall relative to precision scalar → F-β scalar

Let \(y \in \{0, 1\}^n\) be the true labels and \(p \in \{0, 1\}^n\) the predicted labels, let \(\beta > 0\) (the argument b), and let \(P = \operatorname{precisionScore}(y, p)\) and \(R = \operatorname{recallScore}(y, p)\). Then

\[\operatorname{fbeta}(y, p, \beta) = \frac{(1 + \beta^2)\, P\, R}{q + [\, q = 0 \,]}, \qquad q = \beta^2 P + R .\]

The weighted harmonic mean of precision and recall. b = 1 gives \(F_1\); b = 2 counts recall as twice as important; b = 0.5 favours precision.

Intuition \(\beta\) says how much more recall matters than precision. Written in counts, the score weights each miss \(\beta^2\) times as heavily as a false alarm. \(\beta = 1\) is f1; small \(\beta\) moves toward precision and large \(\beta\) toward recall.

Use it when the two errors have known, unequal costs: \(F_2\) for screening, where misses are worse, and \(F_{0.5}\) for recommendations or alerts, where false alarms are worse. Fix b before looking at the results, not after. It is guarded to 0 like f1.

classificationmodel-evaluationthreshold-selection

def fbeta(y, p, b) ≔ (1 + b · b) · precisionScore(y, p) · recallScore(y, p) / safe(b · b · precisionScore(y, p) + recallScore(y, p))

matthews

matthews(y, p) · y true labels vector n · p predicted labels vector n → MCC scalar

\[\operatorname{matthews}(y, p) = \frac{\mathrm{TP} \cdot \mathrm{TN} - \mathrm{FP} \cdot \mathrm{FN}}{\sqrt{u^+ + [\, u^+ = 0 \,]}} .\]

The Matthews correlation coefficient, between −1 and 1: 1 for perfect prediction, 0 for no better than chance, −1 for total disagreement. When any of the four margins is zero the denominator is replaced by 1, which gives 0 in the usual degenerate cases (for instance a classifier that predicts one class only).

Intuition The Pearson correlation between the 0/1 truth and the 0/1 prediction, computed from the four counts. Every cell of the confusion matrix enters it, so it scores high only when the model does well on both classes, and it does not change if you swap which class is called positive.

Use it when classes are imbalanced and you want one honest number for comparing models: 0 is chance, 1 is perfect, and a negative value is worse than chance. It is usually a better single summary than accuracy, and unlike f1 it credits true negatives. The degenerate cases, such as a model that predicts one class only, are guarded to 0.

classificationclass-imbalancemodel-selection

def matthews(y, p)       ≔ matthewsTop(y, p) / sqrt(safe(max(matthewsUnder(y, p), 0)))

Building blocks

matthewsTop(y, p) — Let \(y, p \in \{0, 1\}^n\). Then \(\operatorname{matthewsTop}(y, p) = \mathrm{TP} \cdot \mathrm{TN} - \mathrm{FP} \cdot \mathrm{FN}\). The numerator of the Matthews coefficient.
matthewsUnder(y, p) — Let \(y, p \in \{0, 1\}^n\), \(c = \mathrm{TP} + \mathrm{FP}\) and \(a = \sum_i y_i\). Then \(\operatorname{matthewsUnder}(y, p) = c\, a\, (n - a)(n - c) = (\mathrm{TP} + \mathrm{FP})(\mathrm{TP} + \mathrm{FN})(\mathrm{TN} + \mathrm{FP})(\mathrm{TN} + \mathrm{FN})\). The product of the four margins, under the square root in the denominator.

The area under the curve

The area under the ROC curve equals the probability that a randomly chosen positive scores higher than a randomly chosen negative, with ties counting one half. That is the Mann–Whitney U statistic of the positive class divided by \(n_1 n_2\), so the module reads it from mannWhitney in Prelude.Resampling: no sorting of scores, no threshold sweep and no trapezoid rule.

\[\mathrm{AUC} = \frac{U_{\text{pos}}}{n_{\text{pos}}\, n_{\text{neg}}}\]

auc

auc(g, s) · g one-hot class membership matrix n×2 · s scores vector n → area under the ROC curve scalar

Let \(s \in \mathbb{R}^n\) be scores and let \(g \in \{0, 1\}^{n \times 2}\) be a one-hot membership matrix whose columns, at positions \(0\) and \(1\), mark the negative and the positive class. Let \(r_i\) be the one-based midrank of \(s_i\) among all \(n\) scores, \(n_0 = \sum_i g_{i0}\), \(n_1 = \sum_i g_{i1}\), and let \(U_1 = \sum_i g_{i1}\, r_i - \frac{n_1 (n_1 + 1)}{2}\) be the Mann–Whitney statistic of the positive class. Then

\[\operatorname{auc}(g, s) = \frac{U_1}{n_0\, n_1 + [\, n_0\, n_1 = 0 \,]} .\]

The ROC-AUC of the scores s. The membership matrix has one row per observation; the second column (index 1) marks the positive class and the first marks the negative class. Higher scores should mean “more likely positive”. Tied scores are handled with midranks.

Intuition Pick a random positive and a random negative: AUC is the probability that the positive gets the higher score, with ties counting half. Only the ordering of the scores matters, across every threshold at once, so 0.5 is a random ranking and 1 is perfect separation; calibration and the choice of threshold play no part.

Use it when you compare scoring models before a threshold is chosen, as in credit scoring, fraud or triage. Under extreme imbalance it can look good while the top of the ranking is still poor, so also check precisionScore at your operating threshold. Ranks cost \(O(n^2)\). If only one class is present, the guard returns 0, not 0.5.

classificationrankingmodel-selection

def auc(g, s) ≔ einsum("g,g->", mannWhitney(g, s), positions(mannWhitney(g, s))) / safe(mannWhitneyTotal(g))

The probabilistic scores

These score predicted probabilities rather than hard labels:

\[\text{log loss} = -\frac{1}{n}\sum_i \bigl[\, y_i \ln p_i + (1 - y_i) \ln(1 - p_i) \,\bigr], \qquad \text{Brier} = \frac{1}{n}\sum_i (p_i - y_i)^2\]

logLoss

logLoss(y, p) · y true labels, 0/1 vector n · p predicted probabilities vector n → log loss scalar

Let \(y \in \{0, 1\}^n\) be the true labels and \(p \in [0, 1]^n\) the predicted probabilities of the positive class. Then

\[\operatorname{logLoss}(y, p) = -\frac{1}{n} \sum_{i=1}^{n} \Bigl[\, y_i \ln \max(p_i, 10^{-12}) + (1 - y_i) \ln \max(1 - p_i, 10^{-12}) \,\Bigr] .\]

The binary cross-entropy of the predicted probabilities. Probabilities are floored at \(10^{-12}\) before the logarithm, so a confident wrong prediction gives a large but finite loss.

Intuition The average surprise of the actual outcomes under the predicted probabilities. Predicting 0.9 for something that happens costs about 0.105, and predicting 0.01 costs about 4.6, so confident mistakes are punished hard. It rewards probabilities that are both sharp and honest. Always predicting 0.5 scores \(\ln 2 \approx 0.693\).

Use it when the probabilities themselves drive decisions (pricing risk, bidding, ensembling), or to compare probabilistic classifiers; it is the quantity logistic regression minimises. Unlike auc, it penalises poor calibration. The floor at \(10^{-12}\) caps one certain-and-wrong prediction at about \(27.6/n\). For a bounded score that a few confident errors cannot dominate, use brier.

probability-calibrationclassificationmodel-evaluation

def logLoss(y, p) ≔ logisticLoss(p, y)

brier

brier(y, p) · y true labels, 0/1 vector n · p predicted probabilities vector n → Brier score scalar

Let \(y \in \{0, 1\}^n\) be the true labels and \(p \in [0, 1]^n\) the predicted probabilities. Then

\[\operatorname{brier}(y, p) = \frac{1}{n} \sum_{i=1}^{n} (p_i - y_i)^2 .\]

The mean squared difference between the predicted probability and the outcome. 0 is perfect; always predicting 0.5 scores 0.25.

Intuition Mean squared error applied to probabilities: the squared distance between each forecast and its 0/1 outcome. It rewards forecasts that are both calibrated and sharp, but a confident miss costs at most 1, whereas logLoss can charge far more.

Use it when you evaluate probability forecasts (weather, match outcomes, churn risk) and want a bounded, easily explained score that a few confident errors cannot dominate. 0 is perfect, and a constant forecast of 0.5 scores 0.25. Compare against a forecast that always predicts the base rate before calling a model good.

probability-calibrationforecastingmodel-evaluation

def brier(y, p)   ≔ einsum("i,i->", p − y, p − y) / count(y)

Regression

For real-valued predictions, with residuals \(r_i = y_i - p_i\):

\[\mathrm{MSE} = \frac{1}{n}\sum_i r_i^2, \quad \mathrm{MAE} = \frac{1}{n}\sum_i \lvert r_i \rvert, \quad \mathrm{MAPE} = \frac{1}{n}\sum_i \frac{\lvert r_i \rvert}{\lvert y_i \rvert}, \quad R^2 = 1 - \frac{\sum_i r_i^2}{\sum_i (y_i - \bar{y})^2}\]

\(R^2\) is negative when the predictions are worse than always predicting the mean. That is a correct answer, not an error.

mse

mse(y, p) · y observed values vector n · p predictions vector n → mean squared error scalar

Let \(y \in \mathbb{R}^n\) be the observed values and \(p \in \mathbb{R}^n\) the predictions, with residuals \(r_i = y_i - p_i\). Then

\[\operatorname{mse}(y, p) = \frac{1}{n} \sum_{i=1}^{n} r_i^2 .\]

The mean of the squared residuals.

Intuition The average squared residual. Squaring makes large errors count disproportionately (an error of 10 costs 100 times an error of 1), and the best constant prediction under it is the mean, which is why least-squares models estimate conditional means.

Use it when large errors are especially harmful, or as the training objective and comparison score for least-squares models. Its units are the square of y’s, so report rmse to people. When outliers in y should not dominate the score, use mae.

regressionmodel-evaluation

def mse(y, p)      ≔ einsum("i,i->", residual(y, p), residual(y, p)) / count(y)

rmse

rmse(y, p) · y observed values vector n · p predictions vector n → root mean squared error scalar

Let \(y, p \in \mathbb{R}^n\) be the observed values and the predictions. Then

\[\operatorname{rmse}(y, p) = \sqrt{\max\bigl(\operatorname{mse}(y, p),\, 0\bigr)} = \sqrt{\frac{1}{n} \sum_{i=1}^{n} (y_i - p_i)^2} .\]

The square root of the mean squared error, in the units of y.

Intuition The square root of mse, back in the units of y: roughly the size of a typical error, weighted toward the larger ones. It is never below mae, and a large gap between the two means a few big errors rather than many small ones.

Use it when you report regression or forecast error in natural units (“off by about 3 °C”) and compare models on the same target. It depends on the scale of y, so do not compare it across targets; use mape or r2 for that. The max(…, 0) only protects the square root against rounding.

regressionforecastingmodel-evaluation

def rmse(y, p)     ≔ sqrt(max(mse(y, p), 0))

mae

mae(y, p) · y observed values vector n · p predictions vector n → mean absolute error scalar

Let \(y, p \in \mathbb{R}^n\) be the observed values and the predictions. Then

\[\operatorname{mae}(y, p) = \frac{1}{n} \sum_{i=1}^{n} \lvert y_i - p_i \rvert .\]

The mean of the absolute residuals.

Intuition The average absolute error. Every unit of error costs the same, so one huge miss counts only in proportion to its size, and the best constant prediction under it is the median rather than the mean.

Use it when errors cost linearly (units short in stock, minutes late), the target has outliers you do not want to dominate the score, or you want a figure anyone can read directly. Report it next to rmse: similar values mean errors of even size, a large gap means a few big misses. It is scale dependent, like rmse.

regressionforecastingrobust-statistics

def mae(y, p)      ≔ einsum("i->", abs(residual(y, p))) / count(y)

mape

mape(y, p) · y observed values vector n · p predictions vector n → mean absolute percentage error, as a fraction scalar

Let \(y, p \in \mathbb{R}^n\) be the observed values and the predictions. Then

\[\operatorname{mape}(y, p) = \frac{1}{n} \sum_{i=1}^{n} \frac{\lvert y_i - p_i \rvert}{\max(\lvert y_i \rvert,\, 10^{-12})} .\]

The mean absolute percentage error, returned as a fraction (0.2 means 20%). Observed values of zero are replaced by \(10^{-12}\) in the denominator, so a zero in y produces a very large term rather than a division by zero.

Intuition Each error as a fraction of the actual value, averaged: “off by 8% on average”. Relative errors let you compare accuracy across series of very different size.

Use it when you communicate accuracy to business users or compare forecasts across products of different volume. The result is a fraction (0.08, not 8). It breaks down near zero actuals: a zero in y becomes \(10^{-12}\) in the divisor, so one zero can swamp the average. It is also lopsided: with non-negative forecasts an under-forecast costs at most 100%, while an over-forecast has no limit. For series with zeros, use mae.

forecastingregressionmodel-evaluation

def mape(y, p)     ≔ einsum("i->", abs(residual(y, p)) / max(abs(y), 1e-12)) / count(y)

r2

r2(y, p) · y observed values vector n · p predictions vector n → coefficient of determination scalar

Let \(y, p \in \mathbb{R}^n\) be the observed values and the predictions, let \(\bar{y} = \frac{1}{n} \sum_i y_i\), and let \(T = \sum_{i} (y_i - \bar{y})^2\). Then

\[\operatorname{r2}(y, p) = R^2 = 1 - \frac{\sum_{i=1}^{n} (y_i - p_i)^2}{T + [\, T = 0 \,]} .\]

One minus the residual sum of squares over the total sum of squares about the mean of y. It is at most 1 and may be negative. If y is constant the denominator is replaced by 1.

Intuition The share of y’s variation that the predictions explain, measured against the baseline “always predict the mean of y”. 1 is perfect, 0 is no better than that baseline, and a negative value is worse than it.

Use it when you judge a regression’s fit on a scale-free axis or compare models on the same target. It compares predictions with y directly rather than squaring a correlation, so a biased model is penalised, and on held-out data it can be negative. It says nothing about error size in units, so pair it with rmse. If y is constant, the guard makes it 1 minus the raw sum of squared residuals.

regressionmodel-evaluationmodel-selection

def r2(y, p)        ≔ 1 − einsum("i,i->", residual(y, p), residual(y, p)) / safe(totalSquares(y))

Building blocks

residual(y, p) — Let \(y, p \in \mathbb{R}^n\). Then \(\operatorname{residual}(y, p)_i = y_i - p_i\). The residuals.
totalSquares(y) — Let \(y \in \mathbb{R}^n\) with mean \(\bar{y}\). Then \(\operatorname{totalSquares}(y) = \sum_{i=1}^{n} (y_i - \bar{y})^2\). The total sum of squares of y about its mean; the denominator of r2.

Clustering, judged without labels

The silhouette of a point compares \(a\), its mean distance to its own cluster, with \(b\), its mean distance to the nearest other cluster:

\[s_i = \frac{b_i - a_i}{\max(a_i, b_i)}\]

Values near 1 mean the point sits well inside its cluster; values near 0 mean it lies between two clusters; negative values suggest it belongs elsewhere.

One difference from the textbook definition matters when comparing numbers with other libraries. The mean distance to a point’s own cluster divides by the cluster size \(\lvert C \rvert\) and includes the point’s zero distance to itself, instead of dividing by \(\lvert C \rvert - 1\). This makes \(a\) smaller by the factor \((\lvert C \rvert - 1)/\lvert C \rvert\): half at \(\lvert C \rvert = 2\), under 2% beyond fifty points. Comparing clusterings with different numbers of clusters on the same data is still meaningful, but the values will not match scikit-learn’s exactly, particularly for small clusters. With a single cluster there is no “other” cluster and the result has no meaning.

silhouettes

silhouettes(d, m) · d pairwise distances matrix n×n · m one-hot cluster membership matrix n×K → silhouette of each point vector n

Let \(d \in \mathbb{R}^{n \times n}\) be pairwise distances and let \(m \in \{0, 1\}^{n \times K}\) be a one-hot cluster membership matrix, with cluster sizes \(\lvert C_k \rvert = \sum_j m_{jk}\). Let

\[D_{ik} = \frac{\sum_{j=1}^{n} d_{ij}\, m_{jk}}{\lvert C_k \rvert + [\, \lvert C_k \rvert = 0 \,]}, \qquad a_i = \sum_{k=1}^{K} m_{ik}\, D_{ik}, \qquad b_i = \min_{1 \le k \le K} \bigl(D_{ik} + 10^{300}\, m_{ik}\bigr) .\]

Then

\[\operatorname{silhouettes}(d, m)_i = \frac{b_i - a_i}{\max(a_i, b_i) + [\, \max(a_i, b_i) = 0 \,]}, \qquad i = 1, \dots, n .\]

The silhouette value of every point. d can hold any dissimilarity; m is one-hot, as returned by lowest in Prelude.Kernels for a distance-to-centroid matrix. Every column of m should have at least one member: an empty cluster’s mean distance is guarded to 0, so it would be taken as the nearest other cluster.

Intuition For each point, compare its mean distance to its own cluster (\(a\)) with its mean distance to the closest other cluster (\(b\)). When \(b\) is much larger than \(a\) the point is well placed (near 1); when they are similar it sits on a border (near 0); a negative value means another cluster is closer on average.

Use it when you diagnose a clustering point by point: spotting misassigned points, weak clusters (low average silhouette within one cluster) or outliers. d is a full n×n dissimilarity matrix, so memory grows as \(n^2\). On small clusters the values read higher than scikit-learn’s, because \(a\) includes the point’s own zero distance. It needs at least two non-empty clusters.

clusteringmodel-evaluationanomaly-detection

def silhouettes(d, m)   ≔ (nearestOther(d, m) − ownSpread(d, m)) / safe(max(nearestOther(d, m), ownSpread(d, m)))

silhouette

silhouette(d, m) · d pairwise distances matrix n×n · m one-hot cluster membership matrix n×K → mean silhouette scalar

Let \(d \in \mathbb{R}^{n \times n}\) be pairwise distances and \(m \in \{0, 1\}^{n \times K}\) a one-hot cluster membership matrix, and let \(s = \operatorname{silhouettes}(d, m)\). Then

\[\operatorname{silhouette}(d, m) = \frac{1}{n} \sum_{i=1}^{n} s_i .\]

The mean silhouette over all points, a single score for the whole clustering.

Intuition The mean of the per-point silhouettes: one number for how compact and well separated the clustering is overall. A common rule of thumb reads above about 0.7 as strong structure and below about 0.25 as little real structure.

Use it when you choose the number of clusters (fit several values of \(K\) and keep the highest score) or compare clustering algorithms and distance measures on the same data. It favours compact, convex clusters, so it can undersell good density-based clusterings. Because \(a\) divides by \(\lvert C \rvert\), compare scores within this library rather than against scikit-learn’s.

clusteringmodel-selection

def silhouette(d, m)    ≔ mean(silhouettes(d, m))

Building blocks

clusterSpread(d, m) — Let \(d \in \mathbb{R}^{n \times n}\) and \(m \in \{0, 1\}^{n \times K}\) with \(\lvert C_k \rvert = \sum_j m_{jk}\). Then \(\operatorname{clusterSpread}(d, m)_{ik} = \sum_j d_{ij}\, m_{jk} \,/\, (\lvert C_k \rvert + [\, \lvert C_k \rvert = 0 \,])\). An n×K matrix: each point’s mean distance to each cluster.
ownSpread(d, m) — Let \(d, m\) be as above and \(D = \operatorname{clusterSpread}(d, m)\). Then \(\operatorname{ownSpread}(d, m)_i = a_i = \sum_{k} D_{ik}\, m_{ik}\). Each point’s mean distance to its own cluster.
otherSpreads(d, m) — Let \(d, m\) be as above and \(D = \operatorname{clusterSpread}(d, m)\). Then \(\operatorname{otherSpreads}(d, m)_{ik} = D_{ik} + 10^{300}\, m_{ik}\). The n×K mean distances with the own cluster pushed out of reach.
nearestOther(d, m) — Let \(d, m\) be as above and \(D = \operatorname{clusterSpread}(d, m)\). Then \(\operatorname{nearestOther}(d, m)_i = b_i = \min_{k} \bigl(D_{ik} + 10^{300}\, m_{ik}\bigr)\). Each point’s mean distance to the nearest cluster other than its own.

Examples

Score a binary classifier

Eight observations, three of them positive. The classifier finds two of the three and raises one false alarm.

open Prelude.Metrics

fn accuracyOf(y, p) ↦ accuracy(y, p)
fn precisionOf(y, p) ↦ precisionScore(y, p)
fn recallOf(y, p) ↦ recallScore(y, p)
fn specificityOf(y, p) ↦ specificity(y, p)
fn f1Of(y, p) ↦ f1(y, p)
fn f2Of(y, p) ↦ fbeta(y, p, 2)
fn mccOf(y, p) ↦ matthews(y, p)

using Examples;

var y = new[] { 1d, 1, 1, 0, 0, 0, 0, 0 };
var p = new[] { 1d, 1, 0, 1, 0, 0, 0, 0 };
var acc = Scores.AccuracyOf(y, p);
var prec = Scores.PrecisionOf(y, p);
var rec = Scores.RecallOf(y, p);
var spec = Scores.SpecificityOf(y, p);
var f1 = Scores.F1Of(y, p);
var f2 = Scores.F2Of(y, p);
var mcc = Scores.MccOf(y, p);

Output (verified, rounded to 6 decimal places):

acc = 0.75
prec = 0.666667
rec = 0.666667
spec = 0.8
f1 = 0.666667
f2 = 0.666667
mcc = 0.466667

The counts are TP = 2, FP = 1, FN = 1, TN = 4, so \(\mathrm{MCC} = (8 - 1)/\sqrt{3 \cdot 3 \cdot 5 \cdot 5} = 7/15\). Precision equals recall here, so every \(F_\beta\) equals it too.

A classifier that predicts nothing positive

open Prelude.Metrics

fn precisionOf(y, p) ↦ precisionScore(y, p)

using Examples;

var y = new[] { 1d, 1, 1, 0, 0, 0, 0, 0 };
var none = Scores.PrecisionOf(y, new[] { 0d, 0, 0, 0, 0, 0, 0, 0 });

Output (verified):

prec none = 0

The guarded divisor returns 0 instead of NaN.

ROC-AUC from scores

Three negatives (first column) and four positives (second column).

open Prelude.Metrics

fn aucOf(g, s) ↦ auc(g, s)

using Examples;

var g = new double[,] { {1,0},{1,0},{1,0},{0,1},{0,1},{0,1},{0,1} };
var s = new[] { 0.1d, 0.35, 0.3, 0.8, 0.4, 0.9, 0.2 };
var area = Scores.AucOf(g, s);

Output (verified, rounded to 6 decimal places):

auc = 0.833333

Of the 12 positive–negative pairs, the positive scores higher in 10: the positives 0.8, 0.4 and 0.9 beat all three negatives, and 0.2 beats only 0.1.

Log loss and Brier score

open Prelude.Metrics

fn logLossOf(y, p) ↦ logLoss(y, p)
fn brierOf(y, p) ↦ brier(y, p)

using Examples;

var yb = new[] { 1d, 0, 1, 0 };
var pb = new[] { 0.9d, 0.2, 0.6, 0.4 };
var logLoss = Scores.LogLossOf(yb, pb);
var brier = Scores.BrierOf(yb, pb);

Output (verified, rounded to 6 decimal places):

logloss = 0.337539
brier = 0.0925

The Brier score is \((0.01 + 0.04 + 0.16 + 0.16)/4\).

Regression errors

open Prelude.Metrics

fn mseOf(y, p) ↦ mse(y, p)
fn rmseOf(y, p) ↦ rmse(y, p)
fn maeOf(y, p) ↦ mae(y, p)
fn mapeOf(y, p) ↦ mape(y, p)
fn r2Of(y, p) ↦ r2(y, p)

using Examples;

var yr = new[] { 3d, 5, 2, 7 };
var pr = new[] { 2.5d, 5, 3, 8 };
var mse = Scores.MseOf(yr, pr);
var rmse = Scores.RmseOf(yr, pr);
var mae = Scores.MaeOf(yr, pr);
var mape = Scores.MapeOf(yr, pr);
var r2 = Scores.R2Of(yr, pr);

Output (verified, rounded to 6 decimal places):

mse = 0.5625
rmse = 0.75
mae = 0.625
mape = 0.202381
r2 = 0.847458

The residuals are 0.5, 0, −1, −1; the total sum of squares about the mean 4.25 is 14.75, so \(R^2 = 1 - 2.25/14.75\).

The silhouette of two tight clusters

Four points on a line at 0, 1, 5 and 6 (distances shown as a matrix), clustered as {0, 1} and {5, 6}.

open Prelude.Metrics

fn silhouettesOf(d, m) ↦ silhouettes(d, m)
fn silhouetteOf(d, m) ↦ silhouette(d, m)

using Examples;

var d = new double[,] { {0,1,5,6},{1,0,4,5},{5,4,0,1},{6,5,1,0} };
var m = new double[,] { {1,0},{1,0},{0,1},{0,1} };
var each = Scores.SilhouettesOf(d, m);
var overall = Scores.SilhouetteOf(d, m);

Output (verified, rounded to 6 decimal places):

sils = [0.909091, 0.888889, 0.888889, 0.909091]
sil = 0.89899

For the first point, \(a = (0 + 1)/2 = 0.5\) (its own zero distance included) and \(b = (5 + 6)/2 = 5.5\), so \(s = 5/5.5\). The textbook definition would use \(a = 1\) and give 0.818.