Switch to using Simpson's rule for the integration, and discouple the two

step sizes.

diff --git a/foorank.cpp b/foorank.cpp
index 096c185f3db0ad1ed8159e6f8ac70c3145c0c615..b7d9011194e338df5aed675437f6885144a69e70 100644 (file)
--- a/foorank.cpp
+++ b/foorank.cpp
@@ -5,8 +5,9 @@
 #include <vector>
 #include <algorithm>
 
-// 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<pair<double, double> > 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)));