From 8f9d5fde3ae0b3c659f2a53b57c081476cd7fa14 Mon Sep 17 00:00:00 2001 From: "Steinar H. Gunderson" Date: Thu, 4 Oct 2007 00:07:11 +0200 Subject: [PATCH] Switch to using Simpson's rule for the integration, and discouple the two step sizes. --- foorank.cpp | 40 ++++++++++++++++++++++++++++------------ 1 file changed, 28 insertions(+), 12 deletions(-) diff --git a/foorank.cpp b/foorank.cpp index 096c185..b7d9011 100644 --- a/foorank.cpp +++ b/foorank.cpp @@ -5,8 +5,9 @@ #include #include -// integration step size -static const double step_size = 10.0; +// step sizes +static const double int_step_size = 50.0; +static const double pdf_step_size = 10.0; // rating constant (see below) static const double rating_constant = 455.0; @@ -74,19 +75,34 @@ double prodai(double a) // Set the last parameter to 1.0 if player 1 won, or -1.0 if player 2 won. // In the latter case, ProbScore will be given (r1-r2) instead of (r2-r1). // +static inline double evaluate_int_point(double a, double prodai_precompute, double r1, double mu2, double sigma2, double winfac, double x); + double opponent_rating_pdf(double a, double r1, double mu2, double sigma2, double winfac) { - double sum = 0.0; double prodai_precompute = prodai(a); winfac /= rating_constant; - for (double r2 = 0.0; r2 < 3000.0; r2 += step_size) { - double x = r2 + step_size*0.5; - double probscore = prob_score_real(a, prodai_precompute, (r1 - x)*winfac); - double z = (x - mu2)/sigma2; - double gaussian = exp(-(z*z/2.0)); - sum += step_size * probscore * gaussian; + + int n = int(3000.0 / int_step_size + 0.5); + double h = 3000.0 / double(n); + double sum = evaluate_int_point(a, prodai_precompute, r1, mu2, sigma2, winfac, 0.0); + + for (int i = 1; i < n; i += 2) { + sum += 4.0 * evaluate_int_point(a, prodai_precompute, r1, mu2, sigma2, winfac, i * h); + } + for (int i = 2; i < n; i += 2) { + sum += 2.0 * evaluate_int_point(a, prodai_precompute, r1, mu2, sigma2, winfac, i * h); } - return sum; + sum += evaluate_int_point(a, prodai_precompute, r1, mu2, sigma2, winfac, 3000.0); + + return (h/3.0) * sum; +} + +static inline double evaluate_int_point(double a, double prodai_precompute, double r1, double mu2, double sigma2, double winfac, double x) +{ + double probscore = prob_score_real(a, prodai_precompute, (r1 - x)*winfac); + double z = (x - mu2)/sigma2; + double gaussian = exp(-(z*z/2.0)); + return probscore * gaussian; } // normalize the curve so we know that A ~= 1 @@ -328,13 +344,13 @@ int main(int argc, char **argv) vector > curve; if (score1 == 10) { - for (double r1 = 0.0; r1 < 3000.0; r1 += step_size) { + for (double r1 = 0.0; r1 < 3000.0; r1 += pdf_step_size) { double z = (r1 - mu1) / sigma1; double gaussian = exp(-(z*z/2.0)); curve.push_back(make_pair(r1, gaussian * opponent_rating_pdf(score2, r1, mu2, sigma2, 1.0))); } } else { - for (double r1 = 0.0; r1 < 3000.0; r1 += step_size) { + for (double r1 = 0.0; r1 < 3000.0; r1 += pdf_step_size) { double z = (r1 - mu1) / sigma1; double gaussian = exp(-(z*z/2.0)); curve.push_back(make_pair(r1, gaussian * opponent_rating_pdf(score1, r1, mu2, sigma2, -1.0))); -- 2.39.2