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