X-Git-Url: https://git.sesse.net/?a=blobdiff_plain;f=foosrank.cpp;fp=foosrank.cpp;h=c3bfc88d459c6b35719feda8b3376451cbfd86c0;hb=0797b1180350eae29764a5a282f4ce1bdc69235c;hp=b4f12dc927f2371e03769e2d279eb9fd0505c14c;hpb=6d454fc54830023ce429861a14ef494aa2872445;p=foosball

diff --git a/foosrank.cpp b/foosrank.cpp
index b4f12dc..c3bfc88 100644
--- 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;