+// int(normpdf[mu2, sigma2](t2) * ..., t2=0..3000);
+class OuterIntegralEvaluator {
+private:
+ double theta1, mu2, sigma2, mu_t, sigma_t;
+ int score1, score2;
+ double winfac;
+
+public:
+ OuterIntegralEvaluator(double theta1, double mu2, double sigma2, double mu3, double sigma3, double mu4, double sigma4, int score1, int score2, double winfac)
+ : theta1(theta1), mu2(mu2), sigma2(sigma2), mu_t(mu3 + mu4), sigma_t(sqrt(sigma3*sigma3 + sigma4*sigma4)), score1(score1), score2(score2), winfac(winfac) {}
+
+ double operator() (double theta2) const
+ {
+ double z = (theta2 - mu2) / sigma2;
+ double gaussian = exp(-(z*z/2.0));
+ double r1 = theta1 + theta2;
+ return gaussian * opponent_rating_pdf(score1, score2, r1, mu_t, sigma_t, winfac);
+ }
+};
+
+void compute_new_double_rating(double mu1, double sigma1, double mu2, double sigma2, double mu3, double sigma3, double mu4, double sigma4, int score1, int score2, double &mu, double &sigma)
+{
+ vector<pair<double, double> > curve;
+
+ if (score1 > score2) {
+ 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)));
+ }
+ } 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)));
+ }
+ }
+
+ double mu_est, sigma_est;
+ normalize(curve);
+ estimate_musigma(curve, mu_est, sigma_est);
+ least_squares(curve, mu_est, sigma_est, mu, sigma);
+}
+