From: Steinar H. Gunderson <sesse@debian.org>
Date: Wed, 3 Oct 2007 23:11:06 +0000 (+0200)
Subject: Switch to using Simpson's rule for the moment estimation as well, which seems
X-Git-Url: https://git.sesse.net/?p=foosball;a=commitdiff_plain;h=130255c3b613130f3799dfa17ece66b5c6d62e06

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 b7d9011..f8ecf63 100644
--- 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);