From: Steinar H. Gunderson <sesse@debian.org>
Date: Wed, 3 Oct 2007 21:15:13 +0000 (+0200)
Subject: Fix the initial sigma estimate.
X-Git-Url: https://git.sesse.net/?p=foosball;a=commitdiff_plain;h=7d3cf85d553aad55d17243efcbc76302d3fe9fc8

Fix the initial sigma estimate.
---

diff --git a/foorank.cpp b/foorank.cpp
index 659c611..573a41d 100644
--- a/foorank.cpp
+++ b/foorank.cpp
@@ -215,12 +215,12 @@ void solve3x3(double *A, double *x, double *B)
 }
 
 // Give an OK starting estimate for the least squares, by numerical integration
-// of x*f(x) and x^2 * f(x). Somehow seems to underestimate sigma, though.
+// of statistical moments.
 void estimate_musigma(vector<pair<double, double> > &curve, double &mu_result, double &sigma_result)
 {
-	double mu = 0.0;
-	double sigma = 0.0;
 	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;
@@ -230,17 +230,24 @@ void estimate_musigma(vector<pair<double, double> > &curve, double &mu_result, d
 		double xm = 0.5 * (x0 + x1);
 		double ym = 0.5 * (y0 + y1);
 		sum_area += (x1-x0) * ym;
-		mu += (x1-x0) * xm * ym;
-		sigma += (x1-x0) * xm * xm * ym;
+		ex += (x1-x0) * xm * ym;
+		ex2 += (x1-x0) * xm * xm * ym;
 	}
 
-	mu_result = mu / sum_area;
-	sigma_result = sqrt(sigma) / sum_area;
+	ex /= sum_area;
+	ex2 /= sum_area;
+
+	mu_result = ex;
+	sigma_result = sqrt(ex2 - ex * ex);
 }
 	
 // Find best fit of the data in curves to a Gaussian pdf, based on the
 // given initial estimates. Works by nonlinear least squares, iterating
 // until we're below a certain threshold.
+//
+// Note that the algorithm blows up quite hard if the initial estimate is
+// not good enough. Use estimate_musigma to get a reasonable starting
+// estimate.
 void least_squares(vector<pair<double, double> > &curve, double mu1, double sigma1, double &mu_result, double &sigma_result)
 {
 	double A = 1.0;