Created
January 10, 2020 12:11
-
-
Save rust-play/ad117141e57ad4205dcb4a591ccba823 to your computer and use it in GitHub Desktop.
Code shared from the Rust Playground
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
// EloRating Calculator | |
// http://www.kurims.kyoto-u.ac.jp/~ooura/gamerf.html | |
pub fn gammaln_6(x: f64) -> f64 { | |
// 0 < x < infinity | |
// err = 2.09144255e-18 | |
(x - 0.5) * (x + 6.0975075753906857609437558e+0).ln() - x | |
+ ((((((1.1240582657165407383437999e-2 / (x + 5.0003589884831925541613237e+0) | |
+ 5.0219722703392090725884168e-1) | |
/ (x + 3.9999966300007508932097016e+0) | |
+ 2.0962955353894997733869983e+0) | |
/ (x + 3.0000000467265241458431618e+0) | |
+ 2.2502304753561816836695856e+0) | |
/ (x + 1.9999999996201023058065171e+0) | |
+ 8.5137081316503418312411656e-1) | |
/ (x + 1.0000000000006553243170562e+0) | |
+ 1.2242597732687991784645973e-1) | |
/ x | |
+ 5.6360656189756064967977564e-3) | |
.ln() | |
} | |
pub fn gammaln(x: f64) -> f64 { | |
// to correct gammaln(1.0) == gammaln(2.0) == 0.0 | |
gammaln_6(x) - 1.0 + 1.0 | |
} | |
// http://www.kurims.kyoto-u.ac.jp/~ooura/gamerf.html | |
pub fn erf_6(x: f64) -> f64 { | |
// |x| <= 0.125 | |
// err = 8.5080e-19 | |
let x2 = x * x; | |
(((((-8.492024351869184700200701e-4 * x2 + 5.223878776856181012778436e-3) * &x2 | |
- 2.686616984476423776951305e-2) | |
* x2 | |
+ 1.128379167066213010234749e-1) | |
* x2 | |
- 3.761263890318336015429662e-1) | |
* x2 | |
+ 1.128379167095512573044943e+0) | |
* x | |
} | |
// http://www.kurims.kyoto-u.ac.jp/~ooura/gamerf.html | |
pub fn erfc_8(x: f64) -> f64 { | |
// -infinity < x < infinity | |
// err = 3.63856888e-17 | |
let x2 = x * x; | |
(if x < 6.103997330986881994334338e+0 { | |
2.0 / ((1.269748999651156838985811e+1 * x).exp() + 1.0) | |
} else { | |
0.0 | |
}) + x | |
* (-x2).exp() | |
* (2.963168851992273778336357e-1 / (x2 + 6.121586444955387580549294e-2) | |
+ 1.815811251346370699640955e-1 / (x2 + 5.509427800560020848936831e-1) | |
+ 6.818664514249394930148282e-2 / (x2 + 1.530396620587703969527527e+0) | |
+ 1.569075431619667090378092e-2 / (x2 + 2.999579523113006340465739e+0) | |
+ 2.212901166815175728291522e-3 / (x2 + 4.958677771282467011450533e+0) | |
+ 1.913958130987428643791697e-4 / (x2 + 7.414712510993354068147575e+0) | |
+ 9.710132840105516234434841e-6 / (x2 + 1.047651043565452375901435e+1) | |
+ 1.666424471743077527468010e-7 / (x2 + 1.484555573455979566646185e+1)) | |
} | |
// http://www.kurims.kyoto-u.ac.jp/~ooura/gamerf.html | |
pub fn erfcx_8(x: f64) -> f64 { | |
// -infinity < x < infinity | |
// err = 3.63856888e-17 | |
let x2 = x * x; | |
(if x < 6.103997330986881994334338e+0 { | |
2.0 * x2.exp() / ((1.269748999651156838985811e+1 * x).exp() + 1.0) | |
} else { | |
0.0 | |
}) + x | |
* (2.963168851992273778336357e-1 / (x2 + 6.121586444955387580549294e-2) | |
+ 1.815811251346370699640955e-1 / (x2 + 5.509427800560020848936831e-1) | |
+ 6.818664514249394930148282e-2 / (x2 + 1.530396620587703969527527e+0) | |
+ 1.569075431619667090378092e-2 / (x2 + 2.999579523113006340465739e+0) | |
+ 2.212901166815175728291522e-3 / (x2 + 4.958677771282467011450533e+0) | |
+ 1.913958130987428643791697e-4 / (x2 + 7.414712510993354068147575e+0) | |
+ 9.710132840105516234434841e-6 / (x2 + 1.047651043565452375901435e+1) | |
+ 1.666424471743077527468010e-7 / (x2 + 1.484555573455979566646185e+1)) | |
} | |
pub fn erfcx(x: f64) -> f64 { | |
erfcx_8(x) | |
} | |
pub fn erfc(x: f64) -> f64 { | |
erfc_8(x) | |
} | |
pub fn erf(x: f64) -> f64 { | |
if x.abs() <= 0.125 { | |
erf_6(x) | |
} else if x < 0.0 { | |
erfc(-x) - 1.0 | |
} else { | |
1.0 - erfc(x) | |
} | |
} | |
// https://github.com/jstat/jstat/ | |
// Returns the inverse of the complementary error function | |
pub fn erfcinv(p: f64) -> f64 { | |
if p >= 2.0 { | |
return -100.0; | |
} | |
if p <= 0.0 { | |
return 100.0; | |
} | |
let pp = if p < 1.0 { p } else { 2.0 - p }; | |
let t = (-2.0 * (pp * 0.5).ln()).sqrt(); | |
let mut x = -0.70711 * ((2.30753 + t * 0.27061) / (1.0 + t * (0.99229 + t * 0.04481)) - t); | |
let err = erfc(x) - pp; | |
x += err / (1.12837916709551257 * (-x * x).exp() - x * err); | |
let err = erfc(x) - pp; | |
x += err / (1.12837916709551257 * (-x * x).exp() - x * err); | |
return if p < 1.0 { x } else { -x }; | |
} | |
// https://github.com/jstat/jstat/ | |
// Evaluates the continued fraction for incomplete beta function by modified | |
// Lentz's method. | |
pub fn betacf(x: f64, a: f64, b: f64) -> f64 { | |
let fpmin = 1e-30; | |
let qab = a + b; | |
let qap = a + 1.0; | |
let qam = a - 1.0; | |
let mut c = 1.0; | |
let mut d = 1.0 - qab * x / qap; | |
// These q's will be used in factors that occur in the coefficients | |
if d.abs() < fpmin { | |
d = fpmin; | |
} | |
d = 1.0 / d; | |
let mut h = d; | |
for mi in 1..101 { | |
let m = mi as f64; | |
let m2 = m + m; | |
let aa = m * (b - m) * x / ((qam + m2) * (a + m2)); | |
// One step (the even one) of the recurrence | |
d = 1.0 + aa * d; | |
if d.abs() < fpmin { | |
d = fpmin; | |
} | |
c = 1.0 + aa / c; | |
if c.abs() < fpmin { | |
c = fpmin; | |
} | |
d = 1.0 / d; | |
h *= d * c; | |
let aa = -(a + m) * (qab + m) * x / ((a + m2) * (qap + m2)); | |
// Next step of the recurrence (the odd one) | |
d = 1.0 + aa * d; | |
if d.abs() < fpmin { | |
d = fpmin; | |
} | |
c = 1.0 + aa / c; | |
if c.abs() < fpmin { | |
c = fpmin; | |
} | |
d = 1.0 / d; | |
let del = d * c; | |
h *= del; | |
if (del - 1.0).abs() < 3e-7 { | |
break; | |
} | |
} | |
h | |
} | |
// https://github.com/jstat/jstat/ | |
// Returns the incomplete beta function I_x(a,b) | |
pub fn ibeta(x: f64, a: f64, b: f64) -> f64 { | |
// Factors in front of the continued fraction. | |
let bt = if x == 0.0 || x == 1.0 { | |
0.0 | |
} else { | |
(gammaln(a + b) - gammaln(a) - gammaln(b) + a * x.ln() + b * (1.0 - x).ln()).exp() | |
}; | |
if x < 0.0 || x > 1.0 { | |
std::f64::NAN | |
} else if x < (a + 1.0) / (a + b + 2.0) { | |
// Use continued fraction directly. | |
bt * betacf(x, a, b) / a | |
} else { | |
// else use continued fraction after making the symmetry transformation. | |
1.0 - bt * betacf(1.0 - x, b, a) / b | |
} | |
} | |
// https://github.com/jstat/jstat/ | |
// Returns the inverse of the incomplete beta function | |
pub fn ibetainv(p: f64, a: f64, b: f64) -> f64 { | |
let eps = 1e-8; | |
let a1 = a - 1.0; | |
let b1 = b - 1.0; | |
if p <= 0.0 { | |
return 0.0; | |
} | |
if p >= 1.0 { | |
return 1.0; | |
} | |
let mut x; | |
if a >= 1.0 && b >= 1.0 { | |
let pp = if p < 0.5 { p } else { 1.0 - p }; | |
let t = (-2.0 * pp.ln()).sqrt(); | |
x = (2.30753 + t * 0.27061) / (1.0 + t * (0.99229 + t * 0.04481)) - t; | |
if p < 0.5 { | |
x = -x; | |
} | |
let al = (x * x - 3.0) / 6.0; | |
let h = 2.0 / (1.0 / (2.0 * a - 1.0) + 1.0 / (2.0 * b - 1.0)); | |
let w = (x * (al + h).sqrt() / h) | |
- (1.0 / (2.0 * b - 1.0) - 1.0 / (2.0 * a - 1.0)) * (al + 5.0 / 6.0 - 2.0 / (3.0 * h)); | |
x = a / (a + b * (2.0 * w).exp()); | |
} else { | |
let lna = (a / (a + b)).ln(); | |
let lnb = (b / (a + b)).ln(); | |
let t = (a * lna).exp() / a; | |
let u = (b * lnb).exp() / b; | |
let w = t + u; | |
if p < t / w { | |
x = (a * w * p).powf(1.0 / a); | |
} else { | |
x = 1.0 - (b * w * (1.0 - p)).powf(1.0 / b); | |
} | |
} | |
let afac = -gammaln(a) - gammaln(b) + gammaln(a + b); | |
for j in 0..10 { | |
if x == 0.0 || x == 1.0 { | |
return x; | |
} | |
let err = ibeta(x, a, b) - p; | |
let t = (a1 * x.ln() + b1 * (1.0 - x).ln() + afac).exp(); | |
let u = err / t; | |
let t = u / (1.0 - 0.5 * f64::min(1.0, u * (a1 / x - b1 / (1.0 - x)))); | |
x -= t; | |
if x <= 0.0 { | |
x = 0.5 * (x + t); | |
} | |
if x >= 1.0 { | |
x = 0.5 * (x + t + 1.0); | |
} | |
if t.abs() < eps * x && j > 0 { | |
break; | |
} | |
} | |
x | |
} | |
pub fn elorating(winrate: f64) -> f64 { | |
if winrate <= 0.0 { | |
std::f64::NEG_INFINITY | |
} else if winrate >= 1.0 { | |
std::f64::INFINITY | |
} else { | |
std::f64::consts::LOG10_E * 400.0 * (winrate / (1.0 - winrate)).ln() | |
} | |
} | |
pub fn elorating_inv(rdiff: f64) -> f64 { | |
1.0 / (1.0 + (-std::f64::consts::LN_10 * 0.0025 * rdiff)) | |
} | |
#[allow(dead_code)] | |
pub struct RatingRange { | |
lower_w: f64, | |
lower_l: f64, | |
upper_w: f64, | |
upper_l: f64, | |
lowerrate: f64, | |
upperrate: f64, | |
k: f64, | |
n: f64, | |
alpha: f64, | |
} | |
// https://github.com/tibigame/HandicappedRook/tree/master/tool | |
pub fn clooper_pearson(k: f64, n: f64, alpha: f64) -> RatingRange { | |
let alpha_l = 0.5 - alpha * 0.5; | |
let alpha_h = 0.5 + alpha * 0.5; | |
let nk1 = n - k + 1.0; | |
let k1 = k + 1.0; | |
let nk = n - k; | |
let (lower_w, lower_l, lowerrate): (f64, f64, f64) = if k < nk1 { | |
let lower_w = ibetainv(alpha_l, k, nk1); | |
(lower_w, 1.0 - lower_w, elorating(lower_w)) | |
} else { | |
let lower_l = ibetainv(alpha_h, nk1, k); | |
(1.0 - lower_l, lower_l, -elorating(lower_l)) | |
}; | |
let (upper_w, upper_l, upperrate): (f64, f64, f64) = if k1 < nk { | |
let upper_w = ibetainv(alpha_h, k1, nk); | |
(upper_w, 1.0 - upper_w, elorating(upper_w)) | |
} else { | |
let upper_l = ibetainv(alpha_l, nk, k1); | |
(1.0 - upper_l, upper_l, -elorating(upper_l)) | |
}; | |
RatingRange { | |
lower_w: if lower_w.is_nan() { 0.0 } else { lower_w }, | |
lower_l: if lower_l.is_nan() { 1.0 } else { lower_l }, | |
upper_w: if upper_w.is_nan() { 1.0 } else { upper_w }, | |
upper_l: if upper_l.is_nan() { 0.0 } else { upper_l }, | |
lowerrate: if lowerrate.is_nan() { | |
std::f64::NEG_INFINITY | |
} else { | |
lowerrate | |
}, | |
upperrate: if upperrate.is_nan() { | |
std::f64::INFINITY | |
} else { | |
upperrate | |
}, | |
k: k, | |
n: n, | |
alpha: alpha, | |
} | |
} | |
fn main() { | |
let (win, draw, lose): (f64, f64, f64) = (483.0, 10.0, 384.0); | |
let games = win + draw + lose; | |
let win_hdraw = win + 0.5 * draw; | |
println!("Win-Draw-Lose: {}-{}-{}", win, draw, lose); | |
println!("Games: {}", games); | |
println!("EloRating median: {:+.2}", elorating(win_hdraw / games)); | |
println!("WinRate: {:.2}%", win_hdraw / games * 100.0); | |
println!("DrawRate: {:.2}%", draw / games * 100.0); | |
println!(); | |
println!(" σ ( % ) | EloRating estimated range | Δrange"); | |
println!("---------------------+-----------------------------------------+--------"); | |
let xvec: Vec<f64> = vec![ | |
0.250_000_000, | |
// 0.382_924_922_548 ~= 0.5σ | |
erf(0.5 * std::f64::consts::FRAC_1_SQRT_2), | |
0.500_000_000, | |
// 0.682_689_492_137 ~= 1.0σ | |
erf(1.0 * std::f64::consts::FRAC_1_SQRT_2), | |
0.750_000_000, | |
0.800_000_000, | |
// 0.866_385_597_462 ~= 1.5σ | |
erf(1.5 * std::f64::consts::FRAC_1_SQRT_2), | |
0.900_000_000, | |
0.950_000_000, | |
// 0.954_499_736_104 ~= 2.0σ | |
erf(2.0 * std::f64::consts::FRAC_1_SQRT_2), | |
0.980_000_000, | |
// 0.987_580_669_348 ~= 2.5σ | |
erf(2.5 * std::f64::consts::FRAC_1_SQRT_2), | |
0.990_000_000, | |
// 0.997_300_203_937 ~= 3.0σ | |
erf(3.0 * std::f64::consts::FRAC_1_SQRT_2), | |
0.998_000_000, | |
0.999_000_000, | |
// 0.999_534_741_842 ~= 3.5σ | |
erf(3.5 * std::f64::consts::FRAC_1_SQRT_2), | |
0.999_800_000, | |
0.999_900_000, | |
// 0.999_936_657_516 ~= 4.0σ | |
erf(4.0 * std::f64::consts::FRAC_1_SQRT_2), | |
0.999_980_000, | |
0.999_990_000, | |
// 0.999_993_204_654 ~= 4.5σ | |
erf(4.5 * std::f64::consts::FRAC_1_SQRT_2), | |
0.999_998_000, | |
0.999_999_000, | |
// 0.999_999_426_697 ~= 5.0σ | |
erf(5.0 * std::f64::consts::FRAC_1_SQRT_2), | |
0.999_999_800, | |
0.999_999_900, | |
// 0.999_999_962_021 ~= 5.5σ | |
erf(5.5 * std::f64::consts::FRAC_1_SQRT_2), | |
0.999_999_980, | |
0.999_999_990, | |
0.999_999_998, | |
// 0.999_999_998_027 ~= 6.0σ | |
erf(6.0 * std::f64::consts::FRAC_1_SQRT_2), | |
0.999_999_999, | |
]; | |
for x in xvec { | |
let rating_range = clooper_pearson(win_hdraw, games, x); | |
println!( | |
"{:5.3}σ ({:10.7}%) | {:+8.2} ({:6.2}%) ~ {:+8.2} ({:6.2}%) | {:7.2}", | |
erfcinv(1.0 - x) * std::f64::consts::SQRT_2, | |
100.0 * x, | |
rating_range.lowerrate, | |
rating_range.lower_w * 100.0, | |
rating_range.upperrate, | |
rating_range.upper_w * 100.0, | |
rating_range.upperrate - rating_range.lowerrate, | |
); | |
} | |
} |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment