Switch to using Simpson's rule for the moment estimation as well, which seems

to yield a way better initial estimate.

diff --git a/foorank.cpp b/foorank.cpp
index b7d9011194e338df5aed675437f6885144a69e70..f8ecf63872be9f5abd67af011f531ba595a901b6 100644 (file)
--- a/foorank.cpp
+++ b/foorank.cpp
@@ -237,24 +237,34 @@ void solve3x3(double *A, double *x, double *B)
 // of statistical moments.
 void estimate_musigma(vector<pair<double, double> > &curve, double &mu_result, double &sigma_result)
 {
-       double sum_area = 0.0;
-       double ex = 0.0;
-       double ex2 = 0.0;
-
-       for (unsigned i = 1; i < curve.size(); ++i) {
-               double x1 = curve[i].first;
-               double x0 = curve[i-1].first;
-               double y1 = curve[i].second;
-               double y0 = curve[i-1].second;
-               double xm = 0.5 * (x0 + x1);
-               double ym = 0.5 * (y0 + y1);
-               sum_area += (x1-x0) * ym;
-               ex += (x1-x0) * xm * ym;
-               ex2 += (x1-x0) * xm * xm * ym;
+       double h = (curve.back().first - curve.front().first) / (curve.size() - 1);
+
+       double area = curve.front().second;
+       double ex = curve.front().first * curve.front().second;
+       double ex2 = curve.front().first * curve.front().first * curve.front().second;
+
+       for (unsigned i = 1; i < curve.size() - 1; i += 2) {
+               double x = curve[i].first;
+               double y = curve[i].second;
+               area += 4.0 * y;
+               ex += 4.0 * x * y;
+               ex2 += 4.0 * x * x * y;
        }
+       for (unsigned i = 2; i < curve.size() - 1; i += 2) {
+               double x = curve[i].first;
+               double y = curve[i].second;
+               area += 2.0 * y;
+               ex += 2.0 * x * y;
+               ex2 += 2.0 * x * x * y;
+       }
+       
+       area += curve.back().second;
+       ex += curve.back().first * curve.back().second;
+       ex2 += curve.back().first * curve.back().first * curve.back().second;
 
-       ex /= sum_area;
-       ex2 /= sum_area;
+       area = (h/3.0) * area;
+       ex = (h/3.0) * ex / area;
+       ex2 = (h/3.0) * ex2 / area;
 
        mu_result = ex;
        sigma_result = sqrt(ex2 - ex * ex);