Implement proper probabilities (integrating the entire pdf), so assessments

take into account the standard deviation. This also enables more proper model
evaluation, as it mirrors the "predictive discrepancy" used in the Glicko
paper.

diff --git a/foosrank.cpp b/foosrank.cpp
index d295c212ed04187a5ba42995a2d2a5852bf65150..802464d6e4d81cd8eee49d7ee6c12ddb39bb78c1 100644 (file)
--- a/foosrank.cpp
+++ b/foosrank.cpp
@@ -23,7 +23,6 @@ static const double l_const = M_PI / (2.0 * sqrt(3.0));
 
 using namespace std;
 
-static double prob_score(int k, int a, double rd);
 static double prob_score_real(int k, int a, double binomial, double rd_norm);
 static double prodai(int k, int a);
 static double fac(int x);
@@ -37,6 +36,7 @@ static double sech2(double x)
 }
 #endif
 
+#if 0
 // probability of match ending k-a (k>a) when winnerR - loserR = RD
 //
 //   +inf  
@@ -57,6 +57,7 @@ static double prob_score(int k, int a, double rd)
 {
        return prob_score_real(k, a, prodai(k, a) / fac(k-1), rd/rating_constant);
 }
+#endif
 
 // computes x^a, probably more efficiently than pow(x, a) (but requires that a
 // is n unsigned integer)
@@ -371,7 +372,7 @@ static void least_squares(vector<pair<double, double> > &curve, double mu1, doub
        sigma_result = sigma;
 }
 
-static void compute_new_rating(double mu1, double sigma1, double mu2, double sigma2, int score1, int score2, double &mu, double &sigma)
+static void compute_new_rating(double mu1, double sigma1, double mu2, double sigma2, int score1, int score2, double &mu, double &sigma, double &probability)
 {
        vector<pair<double, double> > curve;
 
@@ -396,13 +397,40 @@ static void compute_new_rating(double mu1, double sigma1, double mu2, double sig
 #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.
+       // (TODO: Use Simpson's rule here.)
+       {
+               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;
+               }
+       
+               // FFT convolution multiplication factor (FFTW computes unnormalized
+               // transforms)
+               sum /= (curve.size() * 2);      
+
+               // pdf normalization factors
+#if USE_LOGISTIC_DISTRIBUTION
+               sum *= M_PI / (sigma1 * 4.0 * sqrt(3.0));
+               sum *= M_PI / (sigma2 * 4.0 * sqrt(3.0));
+#else
+               sum /= (sigma1 * sqrt(2.0 * M_PI));
+               sum /= (sigma2 * sqrt(2.0 * M_PI));
+#endif
+
+               probability = sum;
+       }
+
        double mu_est, sigma_est;
        normalize(curve);
        estimate_musigma(curve, mu_est, sigma_est);
        least_squares(curve, mu_est, sigma_est, mu, sigma);
 }
 
-static 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)
+static 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, double &probability)
 {
        vector<pair<double, double> > curve, newcurve;
        double mu_t = mu3 + mu4;
@@ -455,6 +483,35 @@ static void compute_new_double_rating(double mu1, double sigma1, double mu2, dou
                newcurve.push_back(make_pair(r1, gaussian * sum));
        }
 
+       // Compute the overall probability of the given result, by integrating
+       // 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.)
+       {
+               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;
+               }
+       
+               // FFT convolution multiplication factor (FFTW computes unnormalized
+               // transforms)
+               sum /= (newcurve.size() * 2);   
+
+               // pdf normalization factors
+#if USE_LOGISTIC_DISTRIBUTION
+               sum *= M_PI / (sigma1 * 4.0 * sqrt(3.0));
+               sum *= M_PI / (sigma2 * 4.0 * sqrt(3.0));
+               sum *= M_PI / (sigma_t * 4.0 * sqrt(3.0));
+#else
+               sum /= (sigma1 * sqrt(2.0 * M_PI));
+               sum /= (sigma2 * sqrt(2.0 * M_PI));
+               sum /= (sigma_t * sqrt(2.0 * M_PI));
+#endif
+
+               probability = sum;
+       }
 
        double mu_est, sigma_est;
        normalize(newcurve);
@@ -482,12 +539,12 @@ int main(int argc, char **argv)
                double sigma4 = atof(argv[8]);
                int score1 = atoi(argv[9]);
                int score2 = atoi(argv[10]);
-               double mu, sigma;
-               compute_new_double_rating(mu1, sigma1, mu2, sigma2, mu3, sigma3, mu4, sigma4, score1, score2, mu, sigma);
+               double mu, sigma, probability;
+               compute_new_double_rating(mu1, sigma1, mu2, sigma2, mu3, sigma3, mu4, sigma4, score1, score2, mu, sigma, probability);
                if (score1 > score2) {
-                       printf("%f %f %f\n", mu, sigma, prob_score(score1, score2, mu3+mu4-(mu1+mu2)));
+                       printf("%f %f %f\n", mu, sigma, probability);
                } else {
-                       printf("%f %f %f\n", mu, sigma, prob_score(score2, score1, mu1+mu2-(mu3+mu4)));
+                       printf("%f %f %f\n", mu, sigma, probability);
                }
        } else if (argc > 8) {
                double mu3 = atof(argv[5]);
@@ -495,58 +552,61 @@ 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;
                        double newsigma1_1, newsigma1_2, newsigma2_1, newsigma2_2;
-                       compute_new_double_rating(mu1, sigma1, mu2, sigma2, mu3, sigma3, mu4, sigma4, k, i, newmu1_1, newsigma1_1);
-                       compute_new_double_rating(mu2, sigma2, mu1, sigma1, mu3, sigma3, mu4, sigma4, k, i, newmu1_2, newsigma1_2);
-                       compute_new_double_rating(mu3, sigma3, mu4, sigma4, mu1, sigma1, mu2, sigma2, i, k, newmu2_1, newsigma2_1);
-                       compute_new_double_rating(mu4, sigma4, mu3, sigma3, mu1, sigma1, mu2, sigma2, i, k, newmu2_2, newsigma2_2);
+                       double probability;
+                       compute_new_double_rating(mu1, sigma1, mu2, sigma2, mu3, sigma3, mu4, sigma4, k, i, newmu1_1, newsigma1_1, probability);
+                       compute_new_double_rating(mu2, sigma2, mu1, sigma1, mu3, sigma3, mu4, sigma4, k, i, newmu1_2, newsigma1_2, probability);
+                       compute_new_double_rating(mu3, sigma3, mu4, sigma4, mu1, sigma1, mu2, sigma2, i, k, newmu2_1, newsigma2_1, probability);
+                       compute_new_double_rating(mu4, sigma4, mu3, sigma3, mu1, sigma1, mu2, sigma2, i, k, newmu2_2, newsigma2_2, probability);
                        printf("%u-%u,%f,%+f,%+f,%+f,%+f\n",
-                               k, i, prob_score(k, i, mu3+mu4-(mu1+mu2)), newmu1_1-mu1, newmu1_2-mu2,
+                               k, i, probability, newmu1_1-mu1, newmu1_2-mu2,
                                newmu2_1-mu3, newmu2_2-mu4);
                }
                for (int i = k; i --> 0; ) {
                        double newmu1_1, newmu1_2, newmu2_1, newmu2_2;
                        double newsigma1_1, newsigma1_2, newsigma2_1, newsigma2_2;
-                       compute_new_double_rating(mu1, sigma1, mu2, sigma2, mu3, sigma3, mu4, sigma4, i, k, newmu1_1, newsigma1_1);
-                       compute_new_double_rating(mu2, sigma2, mu1, sigma1, mu3, sigma3, mu4, sigma4, i, k, newmu1_2, newsigma1_2);
-                       compute_new_double_rating(mu3, sigma3, mu4, sigma4, mu1, sigma1, mu2, sigma2, k, i, newmu2_1, newsigma2_1);
-                       compute_new_double_rating(mu4, sigma4, mu3, sigma3, mu1, sigma1, mu2, sigma2, k, i, newmu2_2, newsigma2_2);
+                       double probability;
+                       compute_new_double_rating(mu1, sigma1, mu2, sigma2, mu3, sigma3, mu4, sigma4, k, i, newmu1_1, newsigma1_1, probability);
+                       compute_new_double_rating(mu1, sigma1, mu2, sigma2, mu3, sigma3, mu4, sigma4, i, k, newmu1_1, newsigma1_1, probability);
+                       compute_new_double_rating(mu2, sigma2, mu1, sigma1, mu3, sigma3, mu4, sigma4, i, k, newmu1_2, newsigma1_2, probability);
+                       compute_new_double_rating(mu3, sigma3, mu4, sigma4, mu1, sigma1, mu2, sigma2, k, i, newmu2_1, newsigma2_1, probability);
+                       compute_new_double_rating(mu4, sigma4, mu3, sigma3, mu1, sigma1, mu2, sigma2, k, i, newmu2_2, newsigma2_2, probability);
                        printf("%u-%u,%f,%+f,%+f,%+f,%+f\n",
-                               i, k, prob_score(k, i, mu1+mu2-(mu3+mu4)), newmu1_1-mu1, newmu1_2-mu2,
+                               i, k, probability, newmu1_1-mu1, newmu1_2-mu2,
                                newmu2_1-mu3, newmu2_2-mu4);
                }
        } else if (argc > 6) {
                int score1 = atoi(argv[5]);
                int score2 = atoi(argv[6]);
-               double mu, sigma;
-               compute_new_rating(mu1, sigma1, mu2, sigma2, score1, score2, mu, sigma);
+               double mu, sigma, probability;
+               compute_new_rating(mu1, sigma1, mu2, sigma2, score1, score2, mu, sigma, probability);
 
                if (score1 > score2) {
-                       printf("%f %f %f\n", mu, sigma, prob_score(score1, score2, mu2-mu1));
+                       printf("%f %f %f\n", mu, sigma, probability);
                } else {
-                       printf("%f %f %f\n", mu, sigma, prob_score(score2, score1, mu1-mu2));
+                       printf("%f %f %f\n", mu, sigma, probability);
                }
        } else {
                int k = atoi(argv[5]);
 
                // assess all possible scores
                for (int i = 0; i < k; ++i) {
-                       double newmu1, newmu2, newsigma1, newsigma2;
-                       compute_new_rating(mu1, sigma1, mu2, sigma2, k, i, newmu1, newsigma1);
-                       compute_new_rating(mu2, sigma2, mu1, sigma1, i, k, newmu2, newsigma2);
+                       double newmu1, newmu2, newsigma1, newsigma2, probability;
+                       compute_new_rating(mu1, sigma1, mu2, sigma2, k, i, newmu1, newsigma1, probability);
+                       compute_new_rating(mu2, sigma2, mu1, sigma1, i, k, newmu2, newsigma2, probability);
                        printf("%u-%u,%f,%+f,%+f\n",
-                               k, i, prob_score(k, i, mu2-mu1), newmu1-mu1, newmu2-mu2);
+                               k, i, probability, newmu1-mu1, newmu2-mu2);
                }
                for (int i = k; i --> 0; ) {
-                       double newmu1, newmu2, newsigma1, newsigma2;
-                       compute_new_rating(mu1, sigma1, mu2, sigma2, i, k, newmu1, newsigma1);
-                       compute_new_rating(mu2, sigma2, mu1, sigma1, k, i, newmu2, newsigma2);
+                       double newmu1, newmu2, newsigma1, newsigma2, probability;
+                       compute_new_rating(mu1, sigma1, mu2, sigma2, i, k, newmu1, newsigma1, probability);
+                       compute_new_rating(mu2, sigma2, mu1, sigma1, k, i, newmu2, newsigma2, probability);
                        printf("%u-%u,%f,%+f,%+f\n",
-                               i, k, prob_score(k, i, mu1-mu2), newmu1-mu1, newmu2-mu2);
+                               i, k, probability, newmu1-mu1, newmu2-mu2);
                }
        }