X-Git-Url: https://git.sesse.net/?a=blobdiff_plain;f=foosrank.cpp;h=802464d6e4d81cd8eee49d7ee6c12ddb39bb78c1;hb=8e3849860e66edb4742b3525e1d0fa436cc188b9;hp=d295c212ed04187a5ba42995a2d2a5852bf65150;hpb=f7903bdddb48a163a6902e832247d0f97a8ad664;p=foosball 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); } }