commit de58d01322bf41e731dc76d72d4ff5a6bf055fb8 Author: Anders Olsson Date: Fri Jun 10 15:22:27 2022 +0200 Initial commit. diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..4fffb2f --- /dev/null +++ b/.gitignore @@ -0,0 +1,2 @@ +/target +/Cargo.lock diff --git a/Cargo.toml b/Cargo.toml new file mode 100644 index 0000000..ae16f6d --- /dev/null +++ b/Cargo.toml @@ -0,0 +1,8 @@ +[package] +name = "trueskill-tt" +version = "0.1.0" +edition = "2021" + +# See more keys and their definitions at https://doc.rust-lang.org/cargo/reference/manifest.html + +[dependencies] diff --git a/README.md b/README.md new file mode 100644 index 0000000..18d3515 --- /dev/null +++ b/README.md @@ -0,0 +1,3 @@ +# TrueSkill - Through Time + +Rust port of [TrueSkillThroughTime.py](https://github.com/glandfried/TrueSkillThroughTime.py). diff --git a/src/game.rs b/src/game.rs new file mode 100644 index 0000000..229a2e7 --- /dev/null +++ b/src/game.rs @@ -0,0 +1,540 @@ +use std::cmp::Reverse; +use std::collections::HashSet; + +use crate::{message::DiffMessages, utils, variable::TeamVariable, Gaussian, Player, N00}; + +pub struct Game { + teams: Vec>, + result: Vec, + p_draw: f64, + pub likelihoods: Vec>, + pub evidence: f64, +} + +impl Game { + pub fn new(teams: Vec>, result: Vec, p_draw: f64) -> Self { + if !result.is_empty() { + assert!( + teams.len() == result.len(), + "len(result) and (len(teams) != len(result))" + ); + } + + assert!(p_draw >= 0.0 && p_draw < 1.0, "0.0 <= p_draw < 1.0"); + + if p_draw == 0.0 { + assert!( + result.iter().collect::>().len() == result.len(), + "(p_draw == 0.0) and (len(result) > 0) and (len(set(result)) != len(result))" + ); + } + + let mut this = Self { + teams, + result, + p_draw, + likelihoods: Vec::new(), + evidence: 0.0, + }; + + this.compute_likelihoods(); + + this + } + + fn performance(&self, index: usize) -> Gaussian { + self.teams[index] + .iter() + .fold(N00, |sum, p| sum + p.performance()) + } + + fn partial_evidence(&mut self, d: &[DiffMessages], margin: &[f64], tie: &[bool], e: usize) { + let mu = d[e].prior.mu(); + let sigma = d[e].prior.sigma(); + + if tie[e] { + self.evidence *= utils::cdf(margin[e], mu, sigma) - utils::cdf(-margin[e], mu, sigma) + } else { + self.evidence *= 1.0 - utils::cdf(margin[e], mu, sigma); + } + } + + fn graphical_model( + &mut self, + ) -> ( + Vec, + Vec, + Vec, + Vec, + Vec, + ) { + if self.result.is_empty() { + self.result = (0..self.teams.len() as u16).rev().collect::>(); + } + + let r = &self.result; + let o = sortperm(r); + + let t = (0..self.teams.len()) + .map(|e| TeamVariable { + prior: self.teams[o[e]] + .iter() + .fold(N00, |sum, p| sum + p.performance()), + ..Default::default() + }) + .collect::>(); + + let d = t + .windows(2) + .map(|window| DiffMessages { + prior: window[0].prior - window[1].prior, + ..Default::default() + }) + .collect::>(); + + let tie = (0..d.len()) + .map(|e| r[o[e]] == r[o[e + 1]]) + .collect::>(); + + let margin = (0..d.len()) + .map(|e| { + if self.p_draw == 0.0 { + 0.0 + } else { + let a: f64 = self.teams[o[e]].iter().map(|a| a.beta.powi(2)).sum(); + let b: f64 = self.teams[o[e + 1]].iter().map(|a| a.beta.powi(2)).sum(); + + utils::compute_margin(self.p_draw, (a + b).sqrt()) + } + }) + .collect::>(); + + self.evidence = 1.0; + + (o, t, d, tie, margin) + } + + fn likelihood_analitico(&mut self) -> Vec> { + let (o, t, d, tie, margin) = self.graphical_model(); + + self.partial_evidence(&d, &margin, &tie, 0); + + let d = d[0].prior; + let (mu_trunc, sigma_trunc) = utils::trunc(d.mu(), d.sigma(), margin[0], tie[0]); + + let (delta_div, theta_div_pow2) = if d.sigma() == sigma_trunc { + ( + d.sigma().powi(2) * mu_trunc - sigma_trunc.powi(2) * d.mu(), + f64::INFINITY, + ) + } else { + ( + (d.sigma().powi(2) * mu_trunc - sigma_trunc.powi(2) * d.mu()) + / (d.sigma().powi(2) - sigma_trunc.powi(2)), + (sigma_trunc.powi(2) * d.sigma().powi(2)) + / (d.sigma().powi(2) - sigma_trunc.powi(2)), + ) + }; + + let mut res = Vec::new(); + + for i in 0..t.len() { + let mut team = Vec::new(); + + for j in 0..self.teams[o[i]].len() { + // + let mu = if d.sigma() == sigma_trunc { + 0.0 + } else { + self.teams[o[i]][j].prior.mu() + + (delta_div - d.mu()) * (-1.0f64).powi(if i == 1 { 1 } else { 0 }) + }; + + let sigma_analitico = (theta_div_pow2 + d.sigma().powi(2) + - self.teams[o[i]][j].prior.sigma().powi(2)) + .sqrt(); + + team.push(Gaussian::new(mu, sigma_analitico)); + } + + res.push(team); + } + + if o[0] >= o[1] { + res.swap(0, 1); + } + + res + } + + fn likelihood_teams(&mut self) -> Vec { + let (o, mut t, mut d, tie, margin) = self.graphical_model(); + + let mut step = (f64::INFINITY, f64::INFINITY); + let mut i = 0; + + while ((step.0 > 1e-6) || (step.1 > 1e-6)) && i < 10 { + step = (0.0, 0.0); + + for e in 0..d.len() - 1 { + d[e].prior = t[e].posterior_win() - t[e + 1].posterior_lose(); + + if i == 0 { + let mu = d[e].prior.mu(); + let sigma = d[e].prior.sigma(); + + if tie[e] { + self.evidence *= + utils::cdf(margin[e], mu, sigma) - utils::cdf(-margin[e], mu, sigma) + } else { + self.evidence *= 1.0 - utils::cdf(margin[e], mu, sigma); + } + } + + d[e].likelihood = utils::approx(d[e].prior, margin[e], tie[e]) / d[e].prior; + + let likelihood_lose = t[e].posterior_win() - d[e].likelihood; + let delta = t[e + 1].likelihood_lose.delta(likelihood_lose); + + step = ( + if step.0 > delta.0 { step.0 } else { delta.0 }, + if step.1 > delta.1 { step.1 } else { delta.1 }, + ); + + t[e + 1].likelihood_lose = likelihood_lose; + } + + for e in (1..d.len()).rev() { + d[e].prior = t[e].posterior_win() - t[e + 1].posterior_lose(); + + if i == 0 && e == d.len() - 1 { + self.partial_evidence(&d, &margin, &tie, e); + } + + d[e].likelihood = utils::approx(d[e].prior, margin[e], tie[e]) / d[e].prior; + + let likelihood_win = t[e + 1].posterior_lose() + d[e].likelihood; + let delta = t[e].likelihood_win.delta(likelihood_win); + + step = ( + if step.0 > delta.0 { step.0 } else { delta.0 }, + if step.1 > delta.1 { step.1 } else { delta.1 }, + ); + + t[e].likelihood_win = likelihood_win; + } + + i += 1; + } + + if d.len() == 1 { + self.partial_evidence(&d, &margin, &tie, 0); + + d[0].prior = t[0].posterior_win() - t[1].posterior_lose(); + d[0].likelihood = utils::approx(d[0].prior, margin[0], tie[0]) / d[0].prior; + } + + let t_e = t.len(); + let d_e = d.len(); + + t[0].likelihood_win = t[1].posterior_lose() + d[0].likelihood; + t[t_e - 1].likelihood_lose = t[t_e - 2].posterior_win() - d[d_e - 1].likelihood; + + (0..t.len()) + .map(|e| t[o[e]].likelihood()) + .collect::>() + } + + fn compute_likelihoods(&mut self) { + if self.teams.len() > 2 { + let m_t_ft = self.likelihood_teams(); + + self.likelihoods = (0..self.teams.len()) + .map(|e| { + (0..self.teams[e].len()) + .map(|i| m_t_ft[e] - self.performance(e).exclude(self.teams[e][i].prior)) + .collect::>() + }) + .collect::>(); + } else { + self.likelihoods = self.likelihood_analitico(); + } + } + + pub fn posteriors(&self) -> Vec> { + (0..self.teams.len()) + .map(|e| { + (0..self.teams[e].len()) + .map(|i| self.likelihoods[e][i] * self.teams[e][i].prior) + .collect::>() + }) + .collect::>() + } +} + +fn sortperm(xs: &[u16]) -> Vec { + let mut x = xs.iter().enumerate().collect::>(); + x.sort_unstable_by_key(|(_, x)| Reverse(*x)); + x.into_iter().map(|(i, _)| i).collect() +} + +#[cfg(test)] +mod tests { + use crate::{Gaussian, Player, GAMMA, N_INF}; + + use super::*; + + #[test] + fn test_sortperm() { + assert_eq!(sortperm(&[0, 1, 2, 0]), vec![2, 1, 0, 3]); + } + + #[test] + fn test_1vs1() { + let t_a = Player::new( + Gaussian::new(25.0, 25.0 / 3.0), + 25.0 / 6.0, + 25.0 / 300.0, + N_INF, + ); + + let t_b = Player::new( + Gaussian::new(25.0, 25.0 / 3.0), + 25.0 / 6.0, + 25.0 / 300.0, + N_INF, + ); + + let g = Game::new(vec![vec![t_a], vec![t_b]], vec![0, 1], 0.0); + let p = g.posteriors(); + + let a = p[0][0]; + let b = p[1][0]; + + assert_eq!(a.mu(), 20.79477925612302); + assert_eq!(b.mu(), 29.205220743876975); + assert_eq!(a.sigma(), 7.194481422570443); + + let t_a = Player::new(Gaussian::new(29.0, 1.0), 25.0 / 6.0, GAMMA, N_INF); + let t_b = Player::new(Gaussian::new(25.0, 25.0 / 3.0), 25.0 / 6.0, GAMMA, N_INF); + + let g = Game::new(vec![vec![t_a], vec![t_b]], vec![0, 1], 0.0); + let p = g.posteriors(); + + let a = p[0][0]; + let b = p[1][0]; + + assert_eq!(a.mu(), 28.896475351225412); + assert_eq!(a.sigma(), 0.9966043313004235); + assert_eq!(b.mu(), 32.18921172045737); + assert_eq!(b.sigma(), 6.062063735879715); + + let t_a = Player::new(Gaussian::new(1.139, 0.531), 1.0, 0.2125, N_INF); + let t_b = Player::new(Gaussian::new(15.568, 0.51), 1.0, 0.2125, N_INF); + + let g = Game::new(vec![vec![t_a], vec![t_b]], vec![0, 1], 0.0); + + assert_eq!(g.likelihoods[0][0].sigma(), f64::INFINITY); + assert_eq!(g.likelihoods[1][0].sigma(), f64::INFINITY); + assert_eq!(g.likelihoods[0][0].mu(), 0.0); + } + + #[test] + fn test_1vs1vs1() { + let t_a = Player::new( + Gaussian::new(25.0, 25.0 / 3.0), + 25.0 / 6.0, + 25.0 / 300.0, + N_INF, + ); + let t_b = Player::new( + Gaussian::new(25.0, 25.0 / 3.0), + 25.0 / 6.0, + 25.0 / 300.0, + N_INF, + ); + let t_c = Player::new( + Gaussian::new(25.0, 25.0 / 3.0), + 25.0 / 6.0, + 25.0 / 300.0, + N_INF, + ); + + let g = Game::new(vec![vec![t_a], vec![t_b], vec![t_c]], vec![1, 2, 0], 0.0); + let p = g.posteriors(); + + let a = p[0][0]; + let b = p[1][0]; + let c = p[2][0]; + + assert_eq!(a.mu(), 25.00000000000592); + assert_eq!(a.sigma(), 6.238469796269066); + assert_eq!(b.mu(), 31.31135822129149); + assert_eq!(b.sigma(), 6.69881865477675); + assert_eq!(c.mu(), 18.688641778702593); + assert_eq!(c.sigma(), 6.698818654778007); + + let t_a = Player::new( + Gaussian::new(25.0, 25.0 / 3.0), + 25.0 / 6.0, + 25.0 / 300.0, + N_INF, + ); + let t_b = Player::new( + Gaussian::new(25.0, 25.0 / 3.0), + 25.0 / 6.0, + 25.0 / 300.0, + N_INF, + ); + let t_c = Player::new( + Gaussian::new(25.0, 25.0 / 3.0), + 25.0 / 6.0, + 25.0 / 300.0, + N_INF, + ); + + let g = Game::new(vec![vec![t_a], vec![t_b], vec![t_c]], vec![2, 1, 0], 0.0); + let p = g.posteriors(); + + let a = p[0][0]; + let b = p[1][0]; + let c = p[2][0]; + + assert_eq!(a.mu(), 31.31135822129149); + assert_eq!(a.sigma(), 6.69881865477675); + assert_eq!(b.mu(), 25.00000000000592); + assert_eq!(b.sigma(), 6.238469796269066); + assert_eq!(c.mu(), 18.688641778702593); + assert_eq!(c.sigma(), 6.698818654778007); + + let t_a = Player::new( + Gaussian::new(25.0, 25.0 / 3.0), + 25.0 / 6.0, + 25.0 / 300.0, + N_INF, + ); + let t_b = Player::new( + Gaussian::new(25.0, 25.0 / 3.0), + 25.0 / 6.0, + 25.0 / 300.0, + N_INF, + ); + let t_c = Player::new( + Gaussian::new(25.0, 25.0 / 3.0), + 25.0 / 6.0, + 25.0 / 300.0, + N_INF, + ); + + let g = Game::new(vec![vec![t_a], vec![t_b], vec![t_c]], vec![1, 2, 0], 0.5); + let p = g.posteriors(); + + let a = p[0][0]; + let b = p[1][0]; + let c = p[2][0]; + + assert_eq!(a.mu(), 24.999999999511545); + assert_eq!(a.sigma(), 6.092561128305945); + assert_eq!(b.mu(), 33.37931495595287); + assert_eq!(b.sigma(), 6.483575782278924); + assert_eq!(c.mu(), 16.62068504453558); + assert_eq!(c.sigma(), 6.483575782198122); + } + + #[test] + fn test_1vs1_draw() { + let t_a = Player::new( + Gaussian::new(25.0, 25.0 / 3.0), + 25.0 / 6.0, + 25.0 / 300.0, + N_INF, + ); + + let t_b = Player::new( + Gaussian::new(25.0, 25.0 / 3.0), + 25.0 / 6.0, + 25.0 / 300.0, + N_INF, + ); + + let g = Game::new(vec![vec![t_a], vec![t_b]], vec![0, 0], 0.25); + let p = g.posteriors(); + + let a = p[0][0]; + let b = p[1][0]; + + assert_eq!(a.mu(), 25.0); + assert_eq!(a.sigma(), 6.469480769842277); + assert_eq!(b.mu(), 25.0); + assert_eq!(b.sigma(), 6.469480769842277); + + let t_a = Player::new(Gaussian::new(25.0, 3.0), 25.0 / 6.0, 25.0 / 300.0, N_INF); + + let t_b = Player::new(Gaussian::new(29.0, 2.0), 25.0 / 6.0, 25.0 / 300.0, N_INF); + + let g = Game::new(vec![vec![t_a], vec![t_b]], vec![0, 0], 0.25); + let p = g.posteriors(); + + let a = p[0][0]; + let b = p[1][0]; + + assert_eq!(a.mu(), 25.736001810566616); + assert_eq!(a.sigma(), 2.709956162204711); + assert_eq!(b.mu(), 28.67288808419261); + assert_eq!(b.sigma(), 1.9164711604544398); + } + + #[test] + fn test_1vs1vs1_draw() { + let t_a = Player::new( + Gaussian::new(25.0, 25.0 / 3.0), + 25.0 / 6.0, + 25.0 / 300.0, + N_INF, + ); + let t_b = Player::new( + Gaussian::new(25.0, 25.0 / 3.0), + 25.0 / 6.0, + 25.0 / 300.0, + N_INF, + ); + let t_c = Player::new( + Gaussian::new(25.0, 25.0 / 3.0), + 25.0 / 6.0, + 25.0 / 300.0, + N_INF, + ); + + let g = Game::new(vec![vec![t_a], vec![t_b], vec![t_c]], vec![0, 0, 0], 0.25); + let p = g.posteriors(); + + let a = p[0][0]; + let b = p[1][0]; + let c = p[2][0]; + + assert_eq!(a.mu(), 24.999999999999996); + assert_eq!(a.sigma(), 5.729068664890827); + assert_eq!(b.mu(), 25.000000000000004); + assert_eq!(b.sigma(), 5.707423522433266); + assert_eq!(c.mu(), 24.999999999999996); + assert_eq!(c.sigma(), 5.729068664890825); + + let t_a = Player::new(Gaussian::new(25.0, 3.0), 25.0 / 6.0, 25.0 / 300.0, N_INF); + let t_b = Player::new(Gaussian::new(25.0, 3.0), 25.0 / 6.0, 25.0 / 300.0, N_INF); + let t_c = Player::new(Gaussian::new(29.0, 2.0), 25.0 / 6.0, 25.0 / 300.0, N_INF); + + let g = Game::new(vec![vec![t_a], vec![t_b], vec![t_c]], vec![0, 0, 0], 0.25); + let p = g.posteriors(); + + let a = p[0][0]; + let b = p[1][0]; + let c = p[2][0]; + + assert_eq!(a.mu(), 25.48850755025261); + assert_eq!(a.sigma(), 2.638208444298423); + assert_eq!(b.mu(), 25.51067170990121); + assert_eq!(b.sigma(), 2.6287517663583633); + assert_eq!(c.mu(), 28.555920328820523); + assert_eq!(c.sigma(), 1.8856891308577184); + } +} diff --git a/src/gaussian.rs b/src/gaussian.rs new file mode 100644 index 0000000..6554031 --- /dev/null +++ b/src/gaussian.rs @@ -0,0 +1,220 @@ +use std::fmt; +use std::ops; + +use crate::utils; + +pub const BETA: f64 = 1.0; +pub const MU: f64 = 0.0; +pub const SIGMA: f64 = BETA * 6.0; +pub const GAMMA: f64 = BETA * 0.03; + +pub const N01: Gaussian = Gaussian::new(0.0, 1.0); +pub const N00: Gaussian = Gaussian::new(0.0, 0.0); +pub const N_INF: Gaussian = Gaussian::new(0.0, f64::INFINITY); +pub const N_MS: Gaussian = Gaussian::new(MU, SIGMA); + +#[derive(Clone, Copy, PartialEq, Debug)] +pub struct Gaussian { + mu: f64, + sigma: f64, +} + +impl Gaussian { + pub const fn new(mu: f64, sigma: f64) -> Self { + Gaussian { mu, sigma } + } + + pub fn mu(&self) -> f64 { + self.mu + } + + pub fn sigma(&self) -> f64 { + self.sigma + } + + pub fn tau(&self) -> f64 { + self.mu * self.pi() + } + + pub fn pi(&self) -> f64 { + self.sigma.powi(-2) + } + + pub fn forget(&self, gamma: f64, t: u32) -> Self { + Self::new( + self.mu, + (self.sigma().powi(2) + t as f64 * gamma.powi(2)).sqrt(), + ) + } + + pub fn delta(&self, m: Gaussian) -> (f64, f64) { + ((self.mu() - m.mu()).abs(), (self.sigma() - m.sigma()).abs()) + } + + pub fn exclude(&self, m: Gaussian) -> Gaussian { + Gaussian::new( + self.mu() - m.mu(), + (self.sigma().powi(2) - m.sigma().powi(2)).sqrt(), + ) + } + + /* + def forget(self,gamma,t): + return Gaussian(self.mu, math.sqrt(self.sigma**2 + t*gamma**2)) + + def delta(self, M): + return abs(self.mu - M.mu) , abs(self.sigma - M.sigma) + + def exclude(self, M): + return Gaussian(self.mu - M.mu, math.sqrt(self.sigma**2 - M.sigma**2) ) + + def isapprox(self, M, tol=1e-4): + return (abs(self.mu - M.mu) < tol) and (abs(self.sigma - M.sigma) < tol) + */ +} + +impl Default for Gaussian { + fn default() -> Self { + Gaussian { + mu: MU, + sigma: SIGMA, + } + } +} + +impl fmt::Display for Gaussian { + fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result { + write!(f, "N(mu={:.3}, sigma={:.3})", self.mu, self.sigma) + } +} + +impl ops::Add for Gaussian { + type Output = Gaussian; + + fn add(self, rhs: Gaussian) -> Self::Output { + Gaussian { + mu: self.mu + rhs.mu, + sigma: (self.sigma.powi(2) + rhs.sigma.powi(2)).sqrt(), + } + } +} + +impl ops::Sub for Gaussian { + type Output = Gaussian; + + fn sub(self, rhs: Gaussian) -> Self::Output { + Gaussian { + mu: self.mu - rhs.mu, + sigma: (self.sigma.powi(2) + rhs.sigma.powi(2)).sqrt(), + } + } +} + +impl ops::Mul for Gaussian { + type Output = Gaussian; + + fn mul(self, rhs: Gaussian) -> Self::Output { + let (mu, sigma) = utils::mu_sigma(self.tau() + rhs.tau(), self.pi() + rhs.pi()); + + Gaussian { mu, sigma } + } +} + +impl ops::Div for Gaussian { + type Output = Gaussian; + + fn div(self, rhs: Gaussian) -> Self::Output { + let (mu, sigma) = utils::mu_sigma(self.tau() - rhs.tau(), self.pi() - rhs.pi()); + + Gaussian { mu, sigma } + } +} + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn test_add() { + let n = Gaussian { + mu: 25.0, + sigma: 25.0 / 3.0, + }; + + let m = Gaussian { + mu: 0.0, + sigma: 1.0, + }; + + assert_eq!( + n + m, + Gaussian { + mu: 25.0, + sigma: 8.393118874676116 + } + ); + } + + #[test] + fn test_sub() { + let n = Gaussian { + mu: 25.0, + sigma: 25.0 / 3.0, + }; + + let m = Gaussian { + mu: 1.0, + sigma: 1.0, + }; + + assert_eq!( + n - m, + Gaussian { + mu: 24.0, + sigma: 8.393118874676116 + } + ); + } + + #[test] + fn test_mul() { + let n = Gaussian { + mu: 25.0, + sigma: 25.0 / 3.0, + }; + + let m = Gaussian { + mu: 0.0, + sigma: 1.0, + }; + + assert_eq!( + n * m, + Gaussian { + mu: 0.35488958990536273, + sigma: 0.992876838486922 + } + ); + } + + #[test] + fn test_div() { + let n = Gaussian { + mu: 25.0, + sigma: 25.0 / 3.0, + }; + + let m = Gaussian { + mu: 0.0, + sigma: 1.0, + }; + + assert_eq!( + m / n, + Gaussian { + mu: -0.3652597402597402, + sigma: 1.0072787050317253 + } + ); + } +} diff --git a/src/history.rs b/src/history.rs new file mode 100644 index 0000000..721ac64 --- /dev/null +++ b/src/history.rs @@ -0,0 +1,5 @@ +pub struct History { + mu: f64, + sigma: f64, + gamma: f64, +} diff --git a/src/lib.rs b/src/lib.rs new file mode 100644 index 0000000..c1d1f2c --- /dev/null +++ b/src/lib.rs @@ -0,0 +1,12 @@ +mod game; +mod gaussian; +mod history; +mod message; +mod player; +mod utils; +mod variable; + +pub use game::*; +pub use gaussian::*; +pub use history::*; +pub use player::*; diff --git a/src/main.rs b/src/main.rs new file mode 100644 index 0000000..9c2d3b1 --- /dev/null +++ b/src/main.rs @@ -0,0 +1,9 @@ +use trueskill_tt::*; + +fn main() { + let t_a = Player::new(Gaussian::new(29.0, 1.0), 25.0 / 6.0, GAMMA, N_INF); + let t_b = Player::new(Gaussian::new(25.0, 25.0 / 3.0), 25.0 / 6.0, GAMMA, N_INF); + + let g = Game::new(vec![vec![t_a], vec![t_b]], vec![0, 1], 0.0); + let p = g.posteriors(); +} diff --git a/src/message.rs b/src/message.rs new file mode 100644 index 0000000..e283c9d --- /dev/null +++ b/src/message.rs @@ -0,0 +1,22 @@ +use crate::{Gaussian, N_INF}; + +#[derive(Debug)] +pub struct DiffMessages { + pub prior: Gaussian, + pub likelihood: Gaussian, +} + +impl DiffMessages { + pub fn p(&self) -> Gaussian { + self.prior * self.likelihood + } +} + +impl Default for DiffMessages { + fn default() -> Self { + Self { + prior: N_INF, + likelihood: N_INF, + } + } +} diff --git a/src/player.rs b/src/player.rs new file mode 100644 index 0000000..bb06d89 --- /dev/null +++ b/src/player.rs @@ -0,0 +1,47 @@ +use std::fmt; + +use crate::{Gaussian, BETA, GAMMA, N_INF}; + +#[derive(Debug)] +pub struct Player { + pub prior: Gaussian, + pub beta: f64, + gamma: f64, + prior_draw: Gaussian, +} + +impl Player { + pub fn new(prior: Gaussian, beta: f64, gamma: f64, prior_draw: Gaussian) -> Self { + Player { + prior, + beta, + gamma, + prior_draw, + } + } +} + +impl Player { + pub fn performance(&self) -> Gaussian { + Gaussian::new( + self.prior.mu(), + (self.prior.sigma().powi(2) + self.beta.powi(2)).sqrt(), + ) + } +} + +impl Default for Player { + fn default() -> Self { + Player::new(Gaussian::default(), BETA, GAMMA, N_INF) + } +} + +impl fmt::Display for Player { + fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result { + write!( + f, + "Player({}, beta={:.3}, gamma={:.3})", + self.prior, self.beta, self.gamma + ) + } +} diff --git a/src/utils.rs b/src/utils.rs new file mode 100644 index 0000000..cdf7562 --- /dev/null +++ b/src/utils.rs @@ -0,0 +1,209 @@ +use std::f64::consts::{FRAC_1_SQRT_2, FRAC_2_SQRT_PI, SQRT_2}; + +use crate::Gaussian; + +const SQRT_TAU: f64 = 2.5066282746310002; + +fn erfc(x: f64) -> f64 { + let z = x.abs(); + let t = 1.0 / (1.0 + z / 2.0); + + let a = -0.82215223 + t * 0.17087277; + let b = 1.48851587 + t * a; + let c = -1.13520398 + t * b; + let d = 0.27886807 + t * c; + let e = -0.18628806 + t * d; + let f = 0.09678418 + t * e; + let g = 0.37409196 + t * f; + let h = 1.00002368 + t * g; + + let r = t * (-z * z - 1.26551223 + t * h).exp(); + + if x >= 0.0 { + r + } else { + 2.0 - r + } +} + +fn erfc_inv(mut y: f64) -> f64 { + if y >= 2.0 { + return f64::NEG_INFINITY; + } + + debug_assert!(y >= 0.0, "argument must be nonnegative"); + + if y == 0.0 { + return f64::INFINITY; + } + + if y >= 1.0 { + y = 2.0 - y; + } + + let t = (-2.0 * (y / 2.0).ln()).sqrt(); + + let mut x = FRAC_1_SQRT_2 * ((2.30753 + t * 0.27061) / (1.0 + t * (0.99229 + t * 0.04481)) - t); + + for _ in 0..3 { + let err = erfc(x) - y; + + x += err / (FRAC_2_SQRT_PI * (-(x.powi(2))).exp() - x * err) + } + + if y < 1.0 { + x + } else { + -x + } +} + +fn ppf(p: f64, mu: f64, sigma: f64) -> f64 { + mu - sigma * SQRT_2 * erfc_inv(2.0 * p) +} + +pub(crate) fn mu_sigma(tau: f64, pi: f64) -> (f64, f64) { + if pi > 0.0 { + return (tau / pi, (1.0 / pi).sqrt()); + } + + if pi + 1e-5 < 0.0 { + panic!("sigma should be greater than 0"); + } + + (0.0, f64::INFINITY) +} + +pub(crate) fn cdf(x: f64, mu: f64, sigma: f64) -> f64 { + let z = -(x - mu) / (sigma * SQRT_2); + + 0.5 * erfc(z) +} + +fn pdf(x: f64, mu: f64, sigma: f64) -> f64 { + let normalizer = (SQRT_TAU * sigma).powi(-1); + let functional = (-((x - mu).powi(2)) / (2.0 * sigma.powi(2))).exp(); + + normalizer * functional +} + +/* +def ppf(p, mu, sigma): + return mu - sigma * sqrt2 * erfcinv(2 * p) +*/ + +fn v_w(mu: f64, sigma: f64, margin: f64, tie: bool) -> (f64, f64) { + if !tie { + let alpha = (margin - mu) / sigma; + + let v = pdf(-alpha, 0.0, 1.0) / cdf(-alpha, 0.0, 1.0); + let w = v * (v + (-alpha)); + + (v, w) + } else { + let alpha = (-margin - mu) / sigma; + let beta = (margin - mu) / sigma; + + let v = (pdf(alpha, 0.0, 1.0) - pdf(beta, 0.0, 1.0)) + / (cdf(beta, 0.0, 1.0) - cdf(alpha, 0.0, 1.0)); + let u = (alpha * pdf(alpha, 0.0, 1.0) - beta * pdf(beta, 0.0, 1.0)) + / (cdf(beta, 0.0, 1.0) - cdf(alpha, 0.0, 1.0)); + let w = -(u - v.powi(2)); + + (v, w) + } +} + +pub(crate) fn trunc(mu: f64, sigma: f64, margin: f64, tie: bool) -> (f64, f64) { + let (v, w) = v_w(mu, sigma, margin, tie); + + let mu_trunc = mu + sigma * v; + let sigma_trunc = sigma * (1.0 - w).sqrt(); + + (mu_trunc, sigma_trunc) +} + +pub(crate) fn approx(n: Gaussian, margin: f64, tie: bool) -> Gaussian { + let (mu, sigma) = trunc(n.mu(), n.sigma(), margin, tie); + + Gaussian::new(mu, sigma) +} + +pub(crate) fn compute_margin(p_draw: f64, sd: f64) -> f64 { + ppf(0.5 - p_draw / 2.0, 0.0, sd).abs() +} + +#[cfg(test)] +mod tests { + use crate::{Gaussian, N01}; + + use super::*; + + #[test] + fn test_ppf() { + assert_eq!(ppf(0.3, N01.mu(), N01.sigma()), -0.5244004458961101); + + let n23 = Gaussian::new(2.0, 3.0); + assert_eq!(ppf(0.3, n23.mu(), n23.sigma()), 0.4267986623116695); + } + + #[test] + fn test_cdf() { + assert_eq!(cdf(0.3, N01.mu(), N01.sigma()), 0.6179114097962345); + + let n23 = Gaussian::new(2.0, 3.0); + assert_eq!(cdf(0.3, n23.mu(), n23.sigma()), 0.28547031198297773); + } + + #[test] + fn test_pdf() { + assert_eq!(pdf(0.3, N01.mu(), N01.sigma()), 0.38138781546052414); + + let n23 = Gaussian::new(2.0, 3.0); + assert_eq!(pdf(0.3, n23.mu(), n23.sigma()), 0.11325579143491937); + } + + #[test] + fn test_compute_margin() { + assert_eq!( + compute_margin(0.25, 2.0f64.sqrt() * (25.0 / 6.0)), + 1.8776005988640154 + ); + assert_eq!( + compute_margin(0.25, 3.0f64.sqrt() * (25.0 / 6.0)), + 2.2995817039804787 + ); + assert_eq!( + compute_margin(0.0, 3.0f64.sqrt() * (25.0 / 6.0)), + 2.71348758713328e-7 + ); + assert_eq!( + compute_margin(1.0, 3.0f64.sqrt() * (25.0 / 6.0)), + f64::INFINITY + ); + } + + #[test] + fn test_trunc() { + let g = Gaussian::new(0.0, 1.0); + + assert_eq!( + trunc(g.mu(), g.sigma(), 0.0, false), + (0.7978845368663289, 0.6028103066716792) + ); + + let g = Gaussian::new(0.0, SQRT_2 * (25.0 / 6.0)); + + assert_eq!( + trunc(g.mu(), g.sigma(), 1.8776005988, true), + (0.0, 1.0767055018086311) + ); + + let g = Gaussian::new(12.0, SQRT_2 * (25.0 / 6.0)); + + assert_eq!( + trunc(g.mu(), g.sigma(), 1.8776005988, true), + (0.39009949143595435, 1.034397855300721) + ); + } +} diff --git a/src/variable.rs b/src/variable.rs new file mode 100644 index 0000000..bccee53 --- /dev/null +++ b/src/variable.rs @@ -0,0 +1,38 @@ +use crate::{Gaussian, N_INF}; + +#[derive(Debug)] +pub struct TeamVariable { + pub prior: Gaussian, + pub likelihood_lose: Gaussian, + pub likelihood_win: Gaussian, + pub likelihood_draw: Gaussian, +} + +impl TeamVariable { + pub fn p(&self) -> Gaussian { + self.prior * self.likelihood_lose * self.likelihood_win * self.likelihood_draw + } + + pub fn posterior_win(&self) -> Gaussian { + self.prior * self.likelihood_lose * self.likelihood_draw + } + + pub fn posterior_lose(&self) -> Gaussian { + self.prior * self.likelihood_win * self.likelihood_draw + } + + pub fn likelihood(&self) -> Gaussian { + self.likelihood_win * self.likelihood_lose * self.likelihood_draw + } +} + +impl Default for TeamVariable { + fn default() -> Self { + Self { + prior: N_INF, + likelihood_lose: N_INF, + likelihood_win: N_INF, + likelihood_draw: N_INF, + } + } +}