Skip to content

Instantly share code, notes, and snippets.

@rust-play

rust-play/playground.rs

Created January 10, 2020 12:11
Show Gist options
  • Save rust-play/ad117141e57ad4205dcb4a591ccba823 to your computer and use it in GitHub Desktop.
Save rust-play/ad117141e57ad4205dcb4a591ccba823 to your computer and use it in GitHub Desktop.
Download ZIP
Code shared from the Rust Playground
Raw
playground.rs
// 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