sioul_learn/
platt.rs

1// SPDX-License-Identifier: GPL-3.0-or-later
2// Copyright © 2026 Aurélien Pierre
3
4//! Platt scaling: an SVM's decision values turned into probabilities,
5//! `P(spam | f) = 1 / (1 + exp(A·f + B))`, fitted on scores the SVM gave to
6//! messages it had not learned from (Platt 1999). Fitted by Newton's method
7//! with a backtracking line search, as Lin, Lin & Weng ("A note on Platt's
8//! probabilistic outputs for support vector machines", Machine Learning 68,
9//! 2007) wrote it: their targets, their stopping rule, a stable log-loss.
10
11/// `A` and `B` of the sigmoid.
12#[derive(Debug, Clone, Copy, PartialEq)]
13pub struct Platt {
14    pub a: f64,
15    pub b: f64,
16}
17
18impl Platt {
19    /// The probability that a message with decision value `f` is spam.
20    pub fn probability(&self, f: f64) -> f64 {
21        let t = self.a * f + self.b;
22        // Never exp of a large positive number.
23        if t >= 0.0 { (-t).exp() / (1.0 + (-t).exp()) } else { 1.0 / (1.0 + t.exp()) }
24    }
25
26    /// Whether a higher score always means more likely spam.
27    pub fn is_increasing(&self) -> bool {
28        self.a < 0.0
29    }
30}
31
32/// Fits A and B on `(decision value, is spam)` pairs; None without both kinds.
33pub fn fit(scores: &[(f64, bool)]) -> Option<Platt> {
34    let spam = scores.iter().filter(|(_, s)| *s).count() as f64;
35    let ham = scores.len() as f64 - spam;
36    if spam == 0.0 || ham == 0.0 {
37        return None;
38    }
39    // Targets pulled a little from 0 and 1 (Platt's prior), so separable data stays finite.
40    let high = (spam + 1.0) / (spam + 2.0);
41    let low = 1.0 / (ham + 2.0);
42    let target = |is_spam: bool| if is_spam { high } else { low };
43    let loss = |a: f64, b: f64| -> f64 {
44        scores
45            .iter()
46            .map(|&(f, s)| {
47                let t = f * a + b;
48                let p = target(s);
49                if t >= 0.0 { p * t + (1.0 + (-t).exp()).ln() } else { (p - 1.0) * t + (1.0 + t.exp()).ln() }
50            })
51            .sum()
52    };
53    let (max_iterations, min_step, sigma) = (100, 1e-10, 1e-12);
54    let mut a = 0.0;
55    let mut b = ((ham + 1.0) / (spam + 1.0)).ln();
56    let mut value = loss(a, b);
57    for _ in 0..max_iterations {
58        // Gradient and Hessian (made positive definite by `sigma`).
59        let (mut h11, mut h22, mut h21, mut g1, mut g2) = (sigma, sigma, 0.0, 0.0, 0.0);
60        for &(f, s) in scores {
61            let t = f * a + b;
62            let (p, q) = if t >= 0.0 { ((-t).exp() / (1.0 + (-t).exp()), 1.0 / (1.0 + (-t).exp())) } else { (1.0 / (1.0 + t.exp()), t.exp() / (1.0 + t.exp())) };
63            let d2 = p * q;
64            h11 += f * f * d2;
65            h22 += d2;
66            h21 += f * d2;
67            let d1 = target(s) - p;
68            g1 += f * d1;
69            g2 += d1;
70        }
71        if g1.abs() < 1e-5 && g2.abs() < 1e-5 {
72            break;
73        }
74        let det = h11 * h22 - h21 * h21;
75        let da = -(h22 * g1 - h21 * g2) / det;
76        let db = -(-h21 * g1 + h11 * g2) / det;
77        let gd = g1 * da + g2 * db;
78        let mut step = 1.0;
79        while step >= min_step {
80            let (new_a, new_b) = (a + step * da, b + step * db);
81            let new_value = loss(new_a, new_b);
82            if new_value < value + 0.0001 * step * gd {
83                (a, b, value) = (new_a, new_b, new_value);
84                break;
85            }
86            step /= 2.0;
87        }
88        if step < min_step {
89            // The line search failed: as good as it gets.
90            break;
91        }
92    }
93    Some(Platt { a, b })
94}
95
96#[cfg(test)]
97mod tests {
98    use super::*;
99    use crate::Rng;
100
101    /// Scores drawn from a known sigmoid: the fit finds it again, and it rises.
102    #[test]
103    fn finds_a_known_sigmoid() {
104        let truth = Platt { a: -2.0, b: 0.5 };
105        let mut rng = Rng::new(3);
106        let scores: Vec<(f64, bool)> = (0..20_000)
107            .map(|_| {
108                let f = rng.unit() * 6.0 - 3.0;
109                (f, rng.unit() < truth.probability(f))
110            })
111            .collect();
112        let fitted = fit(&scores).unwrap();
113        assert!((fitted.a - truth.a).abs() < 0.15 && (fitted.b - truth.b).abs() < 0.1, "{fitted:?}");
114        assert!(fitted.is_increasing());
115        let mut last = 0.0;
116        for i in -50..=50 {
117            let p = fitted.probability(f64::from(i) / 10.0);
118            assert!(p > last && (0.0..=1.0).contains(&p));
119            last = p;
120        }
121    }
122
123    /// Fully separated scores stay finite (the targets are pulled from 0 and 1).
124    #[test]
125    fn separated_scores_stay_finite() {
126        let scores: Vec<(f64, bool)> = (0..100).map(|i| (if i % 2 == 0 { 1.0 + f64::from(i) / 100.0 } else { -1.0 - f64::from(i) / 100.0 }, i % 2 == 0)).collect();
127        let fitted = fit(&scores).unwrap();
128        assert!(fitted.a.is_finite() && fitted.b.is_finite() && fitted.is_increasing(), "{fitted:?}");
129        assert!(fitted.probability(1.0) > 0.9 && fitted.probability(-1.0) < 0.1);
130        assert_eq!(fit(&[(1.0, true)]), None, "one kind only");
131    }
132
133    #[test]
134    fn probabilities_never_overflow() {
135        let p = Platt { a: -5.0, b: 0.0 };
136        assert_eq!(p.probability(1e6), 1.0);
137        assert_eq!(p.probability(-1e6), 0.0);
138    }
139}
sioul_learn/platt.rs

sioul_learn/
platt.rs
