git.sesse.net Git - foosball/blob - foorank.cpp

   1 #include <stdio.h>
   2 #include <math.h>
   3 #include <assert.h>
   4
   5 #include <vector>
   6 #include <algorithm>
   7
   8 // step sizes
   9 static const double int_step_size = 50.0;
  10 static const double pdf_step_size = 10.0;
  11
  12 // rating constant (see below)
  13 static const double rating_constant = 455.0;
  14
  15 using namespace std;
  16
  17 double prob_score(double a, double rd);
  18 double prob_score_real(double a, double prodai, double rd_norm);
  19 double prodai(double a);
  20
  21 // probability of match ending 10-a when winnerR - loserR = RD
  22 //
  23 //   +inf
  24 //     /
  25 //    |
  26 //    | Poisson[lambda1, t](a) * Erlang[lambda2, 10](t) dt
  27 //    |
  28 //   /
  29 // -inf
  30 //
  31 // where lambda1 = 1.0, lambda2 = 2^(rd/455)
  32 //
  33 // The constant of 455 is chosen carefully so to match with the
  34 // Glicko/Bradley-Terry assumption that a player rated 400 points over
  35 // his/her opponent will win with a probability of 10/11 =~ 0.90909.
  36 //
  37 double prob_score(double a, double rd)
  38 {
  39         return prob_score_real(a, prodai(a), rd/rating_constant);
  40 }
  41
  42 // Same, but takes in Product(a+i, i=1..9) as an argument in addition to a. Faster
  43 // if you already have that precomputed, and assumes rd is already divided by 455.
  44 double prob_score_real(double a, double prodai, double rd_norm)
  45 {
  46         double nom =
  47                 pow(2.0, -a*rd_norm) * pow(2.0, 10.0*rd_norm) * pow(pow(2.0, -rd_norm) + 1.0, -a)
  48                 * prodai;
  49         double denom = 362880 * pow(1.0 + pow(2.0, rd_norm), 10.0);
  50         return nom/denom;
  51 }
  52
  53 // Calculates Product(a+i, i=1..9) (see above).
  54 double prodai(double a)
  55 {
  56         return (a+1)*(a+2)*(a+3)*(a+4)*(a+5)*(a+6)*(a+7)*(a+8)*(a+9);
  57 }
  58
  59 //
  60 // Computes the integral
  61 //
  62 //   +inf
  63 //    /
  64 //    |
  65 //    | ProbScore[a] (r2-r1) Gaussian[mu2, sigma2] (dr2) dr2
  66 //    |
  67 //   /
  68 // -inf
  69 //
  70 // For practical reasons, -inf and +inf are replaced by 0 and 3000, which
  71 // is reasonable in the this context.
  72 //
  73 // The Gaussian is not normalized.
  74 //
  75 // Set the last parameter to 1.0 if player 1 won, or -1.0 if player 2 won.
  76 // In the latter case, ProbScore will be given (r1-r2) instead of (r2-r1).
  77 //
  78 static inline double evaluate_int_point(double a, double prodai_precompute, double r1, double mu2, double sigma2, double winfac, double x);
  79
  80 double opponent_rating_pdf(double a, double r1, double mu2, double sigma2, double winfac)
  81 {
  82         double prodai_precompute = prodai(a);
  83         winfac /= rating_constant;
  84
  85         int n = int(3000.0 / int_step_size + 0.5);
  86         double h = 3000.0 / double(n);
  87         double sum = evaluate_int_point(a, prodai_precompute, r1, mu2, sigma2, winfac, 0.0);
  88
  89         for (int i = 1; i < n; i += 2) {
  90                 sum += 4.0 * evaluate_int_point(a, prodai_precompute, r1, mu2, sigma2, winfac, i * h);
  91         }
  92         for (int i = 2; i < n; i += 2) {
  93                 sum += 2.0 * evaluate_int_point(a, prodai_precompute, r1, mu2, sigma2, winfac, i * h);
  94         }
  95         sum += evaluate_int_point(a, prodai_precompute, r1, mu2, sigma2, winfac, 3000.0);
  96
  97         return (h/3.0) * sum;
  98 }
  99
 100 static inline double evaluate_int_point(double a, double prodai_precompute, double r1, double mu2, double sigma2, double winfac, double x)
 101 {
 102         double probscore = prob_score_real(a, prodai_precompute, (r1 - x)*winfac);
 103         double z = (x - mu2)/sigma2;
 104         double gaussian = exp(-(z*z/2.0));
 105         return  probscore * gaussian;
 106 }
 107
 108 // normalize the curve so we know that A ~= 1
 109 void normalize(vector<pair<double, double> > &curve)
 110 {
 111         double peak = 0.0;
 112         for (vector<pair<double, double> >::const_iterator i = curve.begin(); i != curve.end(); ++i) {
 113                 peak = max(peak, i->second);
 114         }
 115
 116         double invpeak = 1.0 / peak;
 117         for (vector<pair<double, double> >::iterator i = curve.begin(); i != curve.end(); ++i) {
 118                 i->second *= invpeak;
 119         }
 120 }
 121
 122 // computes matA * matB
 123 void mat_mul(double *matA, unsigned ah, unsigned aw,
 124              double *matB, unsigned bh, unsigned bw,
 125              double *result)
 126 {
 127         assert(aw == bh);
 128         for (unsigned y = 0; y < bw; ++y) {
 129                 for (unsigned x = 0; x < ah; ++x) {
 130                         double sum = 0.0;
 131                         for (unsigned c = 0; c < aw; ++c) {
 132                                 sum += matA[c*ah + x] * matB[y*bh + c];
 133                         }
 134                         result[y*bw + x] = sum;
 135                 }
 136         }
 137 }
 138
 139 // computes matA^T * matB
 140 void mat_mul_trans(double *matA, unsigned ah, unsigned aw,
 141                    double *matB, unsigned bh, unsigned bw,
 142                    double *result)
 143 {
 144         assert(ah == bh);
 145         for (unsigned y = 0; y < bw; ++y) {
 146                 for (unsigned x = 0; x < aw; ++x) {
 147                         double sum = 0.0;
 148                         for (unsigned c = 0; c < ah; ++c) {
 149                                 sum += matA[x*ah + c] * matB[y*bh + c];
 150                         }
 151                         result[y*bw + x] = sum;
 152                 }
 153         }
 154 }
 155
 156 void print3x3(double *M)
 157 {
 158         printf("%f %f %f\n", M[0], M[3], M[6]);
 159         printf("%f %f %f\n", M[1], M[4], M[7]);
 160         printf("%f %f %f\n", M[2], M[5], M[8]);
 161 }
 162
 163 void print3x1(double *M)
 164 {
 165         printf("%f\n", M[0]);
 166         printf("%f\n", M[1]);
 167         printf("%f\n", M[2]);
 168 }
 169
 170 // solves Ax = B by Gauss-Jordan elimination, where A is a 3x3 matrix,
 171 // x is a column vector of length 3 and B is a row vector of length 3.
 172 // Destroys its input in the process.
 173 void solve3x3(double *A, double *x, double *B)
 174 {
 175         // row 1 -= row 0 * (a1/a0)
 176         {
 177                 double f = A[1] / A[0];
 178                 A[1] = 0.0;
 179                 A[4] -= A[3] * f;
 180                 A[7] -= A[6] * f;
 181
 182                 B[1] -= B[0] * f;
 183         }
 184
 185         // row 2 -= row 0 * (a2/a0)
 186         {
 187                 double f = A[2] / A[0];
 188                 A[2] = 0.0;
 189                 A[5] -= A[3] * f;
 190                 A[8] -= A[6] * f;
 191
 192                 B[2] -= B[0] * f;
 193         }
 194
 195         // row 2 -= row 1 * (a5/a4)
 196         {
 197                 double f = A[5] / A[4];
 198                 A[5] = 0.0;
 199                 A[8] -= A[7] * f;
 200
 201                 B[2] -= B[1] * f;
 202         }
 203
 204         // back substitute:
 205
 206         // row 1 -= row 2 * (a7/a8)
 207         {
 208                 double f = A[7] / A[8];
 209                 A[7] = 0.0;
 210
 211                 B[1] -= B[2] * f;
 212         }
 213
 214         // row 0 -= row 2 * (a6/a8)
 215         {
 216                 double f = A[6] / A[8];
 217                 A[6] = 0.0;
 218
 219                 B[0] -= B[2] * f;
 220         }
 221
 222         // row 0 -= row 1 * (a3/a4)
 223         {
 224                 double f = A[3] / A[4];
 225                 A[3] = 0.0;
 226
 227                 B[0] -= B[1] * f;
 228         }
 229
 230         // normalize
 231         x[0] = B[0] / A[0];
 232         x[1] = B[1] / A[4];
 233         x[2] = B[2] / A[8];
 234 }
 235
 236 // Give an OK starting estimate for the least squares, by numerical integration
 237 // of statistical moments.
 238 void estimate_musigma(vector<pair<double, double> > &curve, double &mu_result, double &sigma_result)
 239 {
 240         double h = (curve.back().first - curve.front().first) / (curve.size() - 1);
 241
 242         double area = curve.front().second;
 243         double ex = curve.front().first * curve.front().second;
 244         double ex2 = curve.front().first * curve.front().first * curve.front().second;
 245
 246         for (unsigned i = 1; i < curve.size() - 1; i += 2) {
 247                 double x = curve[i].first;
 248                 double y = curve[i].second;
 249                 area += 4.0 * y;
 250                 ex += 4.0 * x * y;
 251                 ex2 += 4.0 * x * x * y;
 252         }
 253         for (unsigned i = 2; i < curve.size() - 1; i += 2) {
 254                 double x = curve[i].first;
 255                 double y = curve[i].second;
 256                 area += 2.0 * y;
 257                 ex += 2.0 * x * y;
 258                 ex2 += 2.0 * x * x * y;
 259         }
 260
 261         area += curve.back().second;
 262         ex += curve.back().first * curve.back().second;
 263         ex2 += curve.back().first * curve.back().first * curve.back().second;
 264
 265         area = (h/3.0) * area;
 266         ex = (h/3.0) * ex / area;
 267         ex2 = (h/3.0) * ex2 / area;
 268
 269         mu_result = ex;
 270         sigma_result = sqrt(ex2 - ex * ex);
 271 }
 272
 273 // Find best fit of the data in curves to a Gaussian pdf, based on the
 274 // given initial estimates. Works by nonlinear least squares, iterating
 275 // until we're below a certain threshold.
 276 //
 277 // Note that the algorithm blows up quite hard if the initial estimate is
 278 // not good enough. Use estimate_musigma to get a reasonable starting
 279 // estimate.
 280 void least_squares(vector<pair<double, double> > &curve, double mu1, double sigma1, double &mu_result, double &sigma_result)
 281 {
 282         double A = 1.0;
 283         double mu = mu1;
 284         double sigma = sigma1;
 285
 286         // column-major
 287         double matA[curve.size() * 3];  // N x 3
 288         double dbeta[curve.size()];     // N x 1
 289
 290         // A^T * A: 3xN * Nx3 = 3x3
 291         double matATA[3*3];
 292
 293         // A^T * dβ: 3xN * Nx1 = 3x1
 294         double matATdb[3];
 295
 296         double dlambda[3];
 297
 298         for ( ;; ) {
 299                 //printf("A=%f mu=%f sigma=%f\n", A, mu, sigma);
 300
 301                 // fill in A (depends only on x_i, A, mu, sigma -- not y_i)
 302                 for (unsigned i = 0; i < curve.size(); ++i) {
 303                         double x = curve[i].first;
 304
 305                         // df/dA(x_i)
 306                         matA[i + 0 * curve.size()] =
 307                                 exp(-(x-mu)*(x-mu)/(2.0*sigma*sigma));
 308
 309                         // df/dµ(x_i)
 310                         matA[i + 1 * curve.size()] =
 311                                 A * (x-mu)/(sigma*sigma) * matA[i + 0 * curve.size()];
 312
 313                         // df/dσ(x_i)
 314                         matA[i + 2 * curve.size()] =
 315                                 matA[i + 1 * curve.size()] * (x-mu)/sigma;
 316                 }
 317
 318                 // find dβ
 319                 for (unsigned i = 0; i < curve.size(); ++i) {
 320                         double x = curve[i].first;
 321                         double y = curve[i].second;
 322
 323                         dbeta[i] = y - A * exp(- (x-mu)*(x-mu)/(2.0*sigma*sigma));
 324                 }
 325
 326                 // compute a and b
 327                 mat_mul_trans(matA, curve.size(), 3, matA, curve.size(), 3, matATA);
 328                 mat_mul_trans(matA, curve.size(), 3, dbeta, curve.size(), 1, matATdb);
 329
 330                 // solve
 331                 solve3x3(matATA, dlambda, matATdb);
 332
 333                 A += dlambda[0];
 334                 mu += dlambda[1];
 335                 sigma += dlambda[2];
 336
 337                 // terminate when we're down to three digits
 338                 if (fabs(dlambda[0]) <= 1e-3 && fabs(dlambda[1]) <= 1e-3 && fabs(dlambda[2]) <= 1e-3)
 339                         break;
 340         }
 341
 342         mu_result = mu;
 343         sigma_result = sigma;
 344 }
 345
 346 void compute_new_rating(double mu1, double sigma1, double mu2, double sigma2, int score1, int score2, double &mu, double &sigma)
 347 {
 348         vector<pair<double, double> > curve;
 349
 350         if (score1 == 10) {
 351                 for (double r1 = 0.0; r1 < 3000.0; r1 += pdf_step_size) {
 352                         double z = (r1 - mu1) / sigma1;
 353                         double gaussian = exp(-(z*z/2.0));
 354                         curve.push_back(make_pair(r1, gaussian * opponent_rating_pdf(score2, r1, mu2, sigma2, 1.0)));
 355                 }
 356         } else {
 357                 for (double r1 = 0.0; r1 < 3000.0; r1 += pdf_step_size) {
 358                         double z = (r1 - mu1) / sigma1;
 359                         double gaussian = exp(-(z*z/2.0));
 360                         curve.push_back(make_pair(r1, gaussian * opponent_rating_pdf(score1, r1, mu2, sigma2, -1.0)));
 361                 }
 362         }
 363
 364         double mu_est, sigma_est;
 365         normalize(curve);
 366         estimate_musigma(curve, mu_est, sigma_est);
 367         least_squares(curve, mu_est, sigma_est, mu, sigma);
 368 }
 369
 370 int main(int argc, char **argv)
 371 {
 372         double mu1 = atof(argv[1]);
 373         double sigma1 = atof(argv[2]);
 374         double mu2 = atof(argv[3]);
 375         double sigma2 = atof(argv[4]);
 376
 377         if (argc > 5) {
 378                 int score1 = atoi(argv[5]);
 379                 int score2 = atoi(argv[6]);
 380                 double mu, sigma;
 381                 compute_new_rating(mu1, sigma1, mu2, sigma2, score1, score2, mu, sigma);
 382                 printf("%f %f\n", mu, sigma);
 383         } else {
 384                 // assess all possible scores
 385                 for (int i = 0; i <= 9; ++i) {
 386                         double newmu1, newmu2, newsigma1, newsigma2;
 387                         compute_new_rating(mu1, sigma1, mu2, sigma2, 10, i, newmu1, newsigma1);
 388                         compute_new_rating(mu2, sigma2, mu1, sigma1, i, 10, newmu2, newsigma2);
 389                         printf("10-%u,%f,%+f,%+f\n",
 390                                 i, prob_score(i, mu1-mu2), newmu1-mu1, newmu2-mu2);
 391                 }
 392                 for (int i = 10; i --> 0; ) {
 393                         double newmu1, newmu2, newsigma1, newsigma2;
 394                         compute_new_rating(mu1, sigma1, mu2, sigma2, i, 10, newmu1, newsigma1);
 395                         compute_new_rating(mu2, sigma2, mu1, sigma1, 10, i, newmu2, newsigma2);
 396                         printf("%u-10,%f,%+f,%+f\n",
 397                                 i, prob_score(i, mu2-mu1), newmu1-mu1, newmu2-mu2);
 398                 }
 399         }
 400 }
 401