From 6d454fc54830023ce429861a14ef494aa2872445 Mon Sep 17 00:00:00 2001 From: "Steinar H. Gunderson" Date: Sat, 8 Dec 2007 14:43:00 +0100 Subject: [PATCH] Minor integration tweaks. --- foosrank.cpp | 15 +++++++++++---- 1 file changed, 11 insertions(+), 4 deletions(-) diff --git a/foosrank.cpp b/foosrank.cpp index 802464d..b4f12dc 100644 --- a/foosrank.cpp +++ b/foosrank.cpp @@ -400,13 +400,16 @@ static void compute_new_rating(double mu1, double sigma1, double mu2, double sig // 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. - // (TODO: Use Simpson's rule here.) + // (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 += h * h * curve[i].second; + sum += curve[i].second; } + + sum *= h * h; // FFT convolution multiplication factor (FFTW computes unnormalized // transforms) @@ -487,13 +490,17 @@ static void compute_new_double_rating(double mu1, double sigma1, double mu2, dou // the entire resulting pdf. Note that since we're actually evaluating // 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: Use Simpson's rule here.) + // 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 += 4.0 * h * h * h * newcurve[i].second; + sum += newcurve[i].second; } + + sum *= 4.0 * h * h * h; // FFT convolution multiplication factor (FFTW computes unnormalized // transforms) -- 2.39.2