From 8e3849860e66edb4742b3525e1d0fa436cc188b9 Mon Sep 17 00:00:00 2001 From: "Steinar H. Gunderson" Date: Sat, 8 Dec 2007 13:45:53 +0100 Subject: [PATCH] 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. --- foosrank.cpp | 120 ++++++++++++++++++++++++++++++++++++++------------- 1 file changed, 90 insertions(+), 30 deletions(-) diff --git a/foosrank.cpp b/foosrank.cpp index d295c21..802464d 100644 --- 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 > &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 > 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 > 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); } } -- 2.39.2