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 sum_area = 0.0;
 241         double ex = 0.0;
 242         double ex2 = 0.0;
 243
 244         for (unsigned i = 1; i < curve.size(); ++i) {
 245                 double x1 = curve[i].first;
 246                 double x0 = curve[i-1].first;
 247                 double y1 = curve[i].second;
 248                 double y0 = curve[i-1].second;
 249                 double xm = 0.5 * (x0 + x1);
 250                 double ym = 0.5 * (y0 + y1);
 251                 sum_area += (x1-x0) * ym;
 252                 ex += (x1-x0) * xm * ym;
 253                 ex2 += (x1-x0) * xm * xm * ym;
 254         }
 255
 256         ex /= sum_area;
 257         ex2 /= sum_area;
 258
 259         mu_result = ex;
 260         sigma_result = sqrt(ex2 - ex * ex);
 261 }
 262
 263 // Find best fit of the data in curves to a Gaussian pdf, based on the
 264 // given initial estimates. Works by nonlinear least squares, iterating
 265 // until we're below a certain threshold.
 266 //
 267 // Note that the algorithm blows up quite hard if the initial estimate is
 268 // not good enough. Use estimate_musigma to get a reasonable starting
 269 // estimate.
 270 void least_squares(vector<pair<double, double> > &curve, double mu1, double sigma1, double &mu_result, double &sigma_result)
 271 {
 272         double A = 1.0;
 273         double mu = mu1;
 274         double sigma = sigma1;
 275
 276         // column-major
 277         double matA[curve.size() * 3];  // N x 3
 278         double dbeta[curve.size()];     // N x 1
 279
 280         // A^T * A: 3xN * Nx3 = 3x3
 281         double matATA[3*3];
 282
 283         // A^T * dβ: 3xN * Nx1 = 3x1
 284         double matATdb[3];
 285
 286         double dlambda[3];
 287
 288         for ( ;; ) {
 289                 //printf("A=%f mu=%f sigma=%f\n", A, mu, sigma);
 290
 291                 // fill in A (depends only on x_i, A, mu, sigma -- not y_i)
 292                 for (unsigned i = 0; i < curve.size(); ++i) {
 293                         double x = curve[i].first;
 294
 295                         // df/dA(x_i)
 296                         matA[i + 0 * curve.size()] =
 297                                 exp(-(x-mu)*(x-mu)/(2.0*sigma*sigma));
 298
 299                         // df/dµ(x_i)
 300                         matA[i + 1 * curve.size()] =
 301                                 A * (x-mu)/(sigma*sigma) * matA[i + 0 * curve.size()];
 302
 303                         // df/dσ(x_i)
 304                         matA[i + 2 * curve.size()] =
 305                                 matA[i + 1 * curve.size()] * (x-mu)/sigma;
 306                 }
 307
 308                 // find dβ
 309                 for (unsigned i = 0; i < curve.size(); ++i) {
 310                         double x = curve[i].first;
 311                         double y = curve[i].second;
 312
 313                         dbeta[i] = y - A * exp(- (x-mu)*(x-mu)/(2.0*sigma*sigma));
 314                 }
 315
 316                 // compute a and b
 317                 mat_mul_trans(matA, curve.size(), 3, matA, curve.size(), 3, matATA);
 318                 mat_mul_trans(matA, curve.size(), 3, dbeta, curve.size(), 1, matATdb);
 319
 320                 // solve
 321                 solve3x3(matATA, dlambda, matATdb);
 322
 323                 A += dlambda[0];
 324                 mu += dlambda[1];
 325                 sigma += dlambda[2];
 326
 327                 // terminate when we're down to three digits
 328                 if (fabs(dlambda[0]) <= 1e-3 && fabs(dlambda[1]) <= 1e-3 && fabs(dlambda[2]) <= 1e-3)
 329                         break;
 330         }
 331
 332         mu_result = mu;
 333         sigma_result = sigma;
 334 }
 335
 336 int main(int argc, char **argv)
 337 {
 338         double mu1 = atof(argv[1]);
 339         double sigma1 = atof(argv[2]);
 340         double mu2 = atof(argv[3]);
 341         double sigma2 = atof(argv[4]);
 342         int score1 = atoi(argv[5]);
 343         int score2 = atoi(argv[6]);
 344         vector<pair<double, double> > curve;
 345
 346         if (score1 == 10) {
 347                 for (double r1 = 0.0; r1 < 3000.0; r1 += pdf_step_size) {
 348                         double z = (r1 - mu1) / sigma1;
 349                         double gaussian = exp(-(z*z/2.0));
 350                         curve.push_back(make_pair(r1, gaussian * opponent_rating_pdf(score2, r1, mu2, sigma2, 1.0)));
 351                 }
 352         } else {
 353                 for (double r1 = 0.0; r1 < 3000.0; r1 += pdf_step_size) {
 354                         double z = (r1 - mu1) / sigma1;
 355                         double gaussian = exp(-(z*z/2.0));
 356                         curve.push_back(make_pair(r1, gaussian * opponent_rating_pdf(score1, r1, mu2, sigma2, -1.0)));
 357                 }
 358         }
 359
 360         double mu_est, sigma_est, mu, sigma;
 361         normalize(curve);
 362         estimate_musigma(curve, mu_est, sigma_est);
 363         least_squares(curve, mu_est, sigma_est, mu, sigma);
 364         printf("%f %f\n", mu, sigma);
 365 }