Invert the sign in the inner integral, as it's then expressible as a convolution.

diff --git a/foosrank.cpp b/foosrank.cpp
index 6991061fae2973d063d9d879aef7abe70dcf3133..13e9a7914b93539b653a01d55ecdb5acb6019ba8 100644 (file)
--- a/foosrank.cpp
+++ b/foosrank.cpp
@@ -92,7 +92,7 @@ double fac(int x)
 //   +inf
 //    /
 //    |
-//    | ProbScore[a] (r2-r1) Gaussian[mu2, sigma2] (dr2) dr2
+//    | ProbScore[a] (r1-r2) Gaussian[mu2, sigma2] (dr2) dr2
 //    |
 //   /
 // -inf
@@ -116,7 +116,7 @@ public:
                : k(k), a(a), binomial_precompute(binomial_precompute), r1(r1), mu2(mu2), sigma2(sigma2), winfac(winfac) {}
        inline double operator() (double x) const
        {
-               double probscore = prob_score_real(k, a, binomial_precompute, (x - r1)*winfac);
+               double probscore = prob_score_real(k, a, binomial_precompute, (r1 - x)*winfac);
                double z = (x - mu2)/sigma2;
                double gaussian = exp(-(z*z/2.0));
                return probscore * gaussian;
@@ -377,13 +377,13 @@ void compute_new_rating(double mu1, double sigma1, double mu2, double sigma2, in
                for (double r1 = 0.0; r1 < 3000.0; r1 += pdf_step_size) {
                        double z = (r1 - mu1) / sigma1;
                        double gaussian = exp(-(z*z/2.0));
-                       curve.push_back(make_pair(r1, gaussian * opponent_rating_pdf(score1, score2, r1, mu2, sigma2, 1.0)));
+                       curve.push_back(make_pair(r1, gaussian * opponent_rating_pdf(score1, score2, r1, mu2, sigma2, -1.0)));
                }
        } else {
                for (double r1 = 0.0; r1 < 3000.0; r1 += pdf_step_size) {
                        double z = (r1 - mu1) / sigma1;
                        double gaussian = exp(-(z*z/2.0));
-                       curve.push_back(make_pair(r1, gaussian * opponent_rating_pdf(score2, score1, r1, mu2, sigma2, -1.0)));
+                       curve.push_back(make_pair(r1, gaussian * opponent_rating_pdf(score2, score1, r1, mu2, sigma2, 1.0)));
                }
        }
 
@@ -421,13 +421,13 @@ void compute_new_double_rating(double mu1, double sigma1, double mu2, double sig
                for (double r1 = 0.0; r1 < 3000.0; r1 += pdf_step_size) {
                        double z = (r1 - mu1) / sigma1;
                        double gaussian = exp(-(z*z/2.0));
-                       curve.push_back(make_pair(r1, gaussian * simpson_integrate(OuterIntegralEvaluator(r1,mu2,sigma2,mu3,sigma3,mu4,sigma4,score1,score2,0.5), 0.0, 3000.0, int_step_size)));
+                       curve.push_back(make_pair(r1, gaussian * simpson_integrate(OuterIntegralEvaluator(r1,mu2,sigma2,mu3,sigma3,mu4,sigma4,score1,score2,-0.5), 0.0, 3000.0, int_step_size)));
                }
        } else {
                for (double r1 = 0.0; r1 < 3000.0; r1 += pdf_step_size) {
                        double z = (r1 - mu1) / sigma1;
                        double gaussian = exp(-(z*z/2.0));
-                       curve.push_back(make_pair(r1, gaussian * simpson_integrate(OuterIntegralEvaluator(r1,mu2,sigma2,mu3,sigma3,mu4,sigma4,score2,score1,-0.5), 0.0, 3000.0, int_step_size)));
+                       curve.push_back(make_pair(r1, gaussian * simpson_integrate(OuterIntegralEvaluator(r1,mu2,sigma2,mu3,sigma3,mu4,sigma4,score2,score1,0.5), 0.0, 3000.0, int_step_size)));
                }
        }