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