Now Simpson's rule works like a charm... Odd.

diff --git a/foosrank.cpp b/foosrank.cpp
index b4f12dc927f2371e03769e2d279eb9fd0505c14c..c3bfc88d459c6b35719feda8b3376451cbfd86c0 100644 (file)
--- a/foosrank.cpp
+++ b/foosrank.cpp
@@ -396,20 +396,21 @@ static void compute_new_rating(double mu1, double sigma1, double mu2, double sig
                curve[i].second *= gaussian;
 #endif
        }
-
+       
        // Compute the overall probability of the given result, by integrating
        // the entire resulting pdf. Note that since we're actually evaluating
        // a double integral, we'll need to multiply by h² instead of h.
-       // (For some reason, Simpson's rule gives markedly worse accuracy here.
-       // Probably related to h² somehow?)
        {
                double h = (curve.back().first - curve.front().first) / (curve.size() - 1);
-               double sum = 0.0;
-               for (unsigned i = 0; i < curve.size(); ++i) {
-                       sum += curve[i].second;
+               double sum = curve.front().second;
+               for (unsigned i = 1; i < curve.size() - 1; i += 2) {
+                       sum += 4.0 * curve[i].second;
                }
-
-               sum *= h * h;
+               for (unsigned i = 2; i < curve.size() - 1; i += 2) {
+                       sum += 2.0 * curve[i].second;
+               }
+               sum += curve.back().second;
+               sum *= h * h / 3.0;
        
                // FFT convolution multiplication factor (FFTW computes unnormalized
                // transforms)
@@ -491,16 +492,18 @@ static void compute_new_double_rating(double mu1, double sigma1, double mu2, dou
        // a triple integral, we'll need to multiply by 4h³ (no idea where the
        // 4 factor comes from, probably from the 0..6000 range somehow) instead
        // of h.
-       // (TODO: Evaluate Simpson's rule here, although it's probably even worse
-       // than for the single case.)
        {
                double h = (newcurve.back().first - newcurve.front().first) / (newcurve.size() - 1);
-               double sum = 0.0;
-               for (unsigned i = 0; i < newcurve.size(); ++i) {
-                       sum += newcurve[i].second;
+               double sum = newcurve.front().second;
+               for (unsigned i = 1; i < newcurve.size() - 1; i += 2) {
+                       sum += 4.0 * newcurve[i].second;
+               }
+               for (unsigned i = 2; i < newcurve.size() - 1; i += 2) {
+                       sum += 2.0 * newcurve[i].second;
                }
+               sum += newcurve.back().second;
 
-               sum *= 4.0 * h * h * h;
+               sum *= 4.0 * h * h * h / 3.0;
        
                // FFT convolution multiplication factor (FFTW computes unnormalized
                // transforms)
@@ -559,7 +562,7 @@ int main(int argc, char **argv)
                double mu4 = atof(argv[7]);
                double sigma4 = atof(argv[8]);
                int k = atoi(argv[9]);
-       
+
                // assess all possible scores
                for (int i = 0; i < k; ++i) {
                        double newmu1_1, newmu1_2, newmu2_1, newmu2_2;