git.sesse.net Git - wloh/blob - bayeswf.cpp

   1 #include <stdio.h>
   2 #include <math.h>
   3 #include <string.h>
   4 #include <stdlib.h>
   5
   6 #include <map>
   7 #include <vector>
   8 #include <string>
   9 #include <algorithm>
  10
  11 using namespace std;
  12
  13 #define PRIOR_MU 1500
  14 #define PRIOR_WEIGHT 1.0
  15 #define MAX_PLAYERS 4096
  16 #define DUMP_RAW 0
  17
  18 float mu[MAX_PLAYERS];
  19 float sigma[MAX_PLAYERS];
  20 float prior_sigma = 70.0f;
  21
  22 #define EPSILON 1e-3
  23
  24 /*
  25  * L(mu_vec, sigma_vec, matches) = product[ L(mu_A, sigma_A, mu_B, sigma_B, score_AB - score_BA) ]
  26  * log-likelihood = sum[ log( L(mu_A, sigma_A, mu_B, sigma_B, score_AB - score_BA) ) ]
  27  *
  28  * L(mu1, sigma1, mu2, sigma2, score2 - score1) = sigmoid(mu2 - mu1, sqrt(sigma1² + sigma2²), (score2 - score1))
  29  *
  30  * pdf := 1/(sigma * sqrt(2*Pi)) * exp(-(x - mu)^2 / (2 * sigma^2));
  31  * pdfs := subs({ mu = mu1 - mu2, sigma = sqrt(sigma1^2 + sigma2^2) }, pdf);
  32  * diff(log(pdfs), mu1);
  33  */
  34
  35 struct match {
  36         int other_player;
  37         int margin;
  38         float weight;
  39 };
  40 map<int, vector<match> > matches_for_player;
  41
  42 void dump_scores(const vector<string> &players, const float *mu, const float *sigma, int num_players)
  43 {
  44 #if 0
  45         for (int i = 0; i < num_players; ++i) {
  46                 printf("%s=[%5.1f, %4.1f] ", players[i].c_str(), mu[i], sigma[i]);
  47         }
  48         printf("\n");
  49 #elif 0
  50         for (int i = 0; i < num_players; ++i) {
  51                 printf("%5.1f ", mu[i]);
  52         }
  53         printf("\n");
  54 #else
  55         for (int i = 0; i < num_players; ++i) {
  56                 printf("%f %s\n", mu[i], players[i].c_str());
  57         }
  58 #endif
  59 }
  60
  61 /*
  62  * diff(logL, mu1) = -w * (mu1 - mu2 - x) / sigma_c^2
  63  * maximizer for mu1 is given by: sum_i[ (w_i/sigma_c_i)^2 (mu1 - mu2_i - x_i) ] = 0
  64  *                                sum_i[ (w_i/sigma_c_i)^2 mu1 ] = sum_i [ (w_i/sigma_c_i)^2 ( mu2_i + x_i ) ]
  65  *                                mu1 = sum_i [ (w_i/sigma_c_i)^2 ( mu2_i + x_i ) ] / sum_i[ (w_i/sigma_c_i)^2 ]
  66  */
  67 void update_mu(float *mu, float *sigma, int player_num, const vector<match> &matches)
  68 {
  69         if (matches.empty()) {
  70                 return;
  71         }
  72
  73         float nom = 0.0f, denom = 0.0f;
  74
  75         // Prior.
  76         {
  77                 float inv_sigma2 = 1.0f / (prior_sigma * prior_sigma);
  78                 nom += PRIOR_WEIGHT * PRIOR_MU * inv_sigma2;
  79                 denom += PRIOR_WEIGHT * inv_sigma2;
  80         }
  81
  82         // All matches.
  83         for (unsigned i = 0; i < matches.size(); ++i) {
  84                 float sigma1 = sigma[player_num];
  85                 float sigma2 = sigma[matches[i].other_player];
  86                 float inv_sigma_c2 = matches[i].weight / (sigma1 * sigma1 + sigma2 * sigma2);
  87                 float x = matches[i].margin; // / 70.0f;
  88
  89                 nom += (mu[matches[i].other_player] + x) * inv_sigma_c2;
  90                 denom += inv_sigma_c2;
  91         }
  92         mu[player_num] = nom / denom;
  93 }
  94
  95 void dump_raw(const float *mu, const float *sigma, int num_players)
  96 {
  97         for (int i = 0; i < num_players; ++i) {
  98                 for (unsigned j = 0; j < matches_for_player[i].size(); ++j) {
  99                         const match& m = matches_for_player[i][j];
 100
 101                         // Only count each match once.
 102                         if (m.other_player <= i) {
 103                                 continue;
 104                         }
 105
 106                         float mu1 = mu[i];
 107                         float mu2 = mu[m.other_player];
 108                         float sigma1 = sigma[i];
 109                         float sigma2 = sigma[m.other_player];
 110                         float sigma = sqrt(sigma1 * sigma1 + sigma2 * sigma2);
 111                         float mu = mu1 - mu2;
 112                         float x = m.margin;
 113                         float w = m.weight;
 114
 115                         printf("%f %f\n", (x - mu) / sigma, w);
 116                 }
 117         }
 118 }
 119
 120 /*
 121  * diff(logL, sigma1) = sigma1 (-sigma1² - sigma2² + (x - mu)²) / sigma_c²
 122  * maximizer for sigma1 is given by: sum_i[ (1/sigma_c_i)² sigma1 ((x - mu)² - (sigma1² + sigma2²) ] = 0
 123  *                                   sum_i[ (x - mu)² - sigma1² - sigma2² ] = 0                                  |: sigma1 != 0, sigma2 != 0
 124  *                                   sum_i[ (x - mu)² - sigma2² ] = sum[ sigma1² ]
 125  *                                   sigma1 = sqrt( sum_i[ (x - mu)² - sigma2² ] / N )
 126  */
 127 void update_sigma(float *mu, float *sigma, int player_num, const vector<match> &matches)
 128 {
 129         if (matches.size() < 2) {
 130                 return;
 131         }
 132
 133         float sum = 0.0f;
 134         for (unsigned i = 0; i < matches.size(); ++i) {
 135                 float mu1 = mu[player_num];
 136                 float mu2 = mu[matches[i].other_player];
 137                 float mu = mu1 - mu2;
 138                 float sigma2 = sigma[matches[i].other_player];
 139                 float x = matches[i].margin;
 140
 141                 //fprintf(stderr, "x=%f mu=%f sigma2=%f   add %f-%f = %f\n", x, mu, sigma2, (x-mu)*(x-mu), sigma2*sigma2, (x - mu) * (x - mu) - sigma2 * sigma2);
 142                 sum += (x - mu) * (x - mu) - sigma2 * sigma2;
 143         }
 144
 145         if (sum <= 0) {
 146                 return;
 147         }
 148         //fprintf(stderr, "sum=%f\n", sum);
 149         sigma[player_num] = sqrt(sum / matches.size());
 150 }
 151
 152 /*
 153  * diff(logL, sigma) = w ( (x - mu)² - sigma² ) / sigma³
 154  * maximizer for sigma is given by: sum_i[ (w_i/sigma)³ ((x - mu)² - sigma²) ] = 0
 155  *                                   sum_i[ w_i ( (x - mu)² - sigma² ) ] = 0                            |: sigma != 0
 156  *                                   sum_i[ w_i (x - mu)² ] = sum[ w_i sigma² ]
 157  *                                   sigma = sqrt( sum_i[ w_i (x - mu)² ] / sum[w_i] )
 158  */
 159 void update_global_sigma(float *mu, float *sigma, int num_players)
 160 {
 161         float nom = 0.0f, denom = 0.0f;
 162         for (int i = 0; i < num_players; ++i) {
 163                 for (unsigned j = 0; j < matches_for_player[i].size(); ++j) {
 164                         const match& m = matches_for_player[i][j];
 165
 166                         // Only count each match once.
 167                         if (m.other_player <= i) {
 168                                 continue;
 169                         }
 170
 171                         float mu1 = mu[i];
 172                         float mu2 = mu[m.other_player];
 173                         float mu = mu1 - mu2;
 174                         float x = m.margin;
 175                         float w = m.weight;
 176
 177                         nom += w * ((x - mu) * (x - mu));
 178                         denom += w;
 179                 }
 180         }
 181
 182         float best_sigma = sqrt(nom / denom) / sqrt(2.0f);  // Divide evenly between the two players.
 183         for (int i = 0; i < num_players; ++i) {
 184                 sigma[i] = best_sigma;
 185         }
 186 }
 187
 188 /*
 189  * diff(priorlogL, sigma) = w ( (x - mu)² - sigma² ) / sigma³
 190  * maximizer for sigma is given by: sum_i[ (w_i/sigma)³ ((x - mu)² - sigma²) ] = 0
 191  *                                   sum_i[ w_i ( (x - mu)² - sigma² ) ] = 0                            |: sigma != 0
 192  *                                   sum_i[ w_i (x - mu)² ] = sum[ w_i sigma² ]
 193  *                                   sigma = sqrt( sum_i[ w_i (x - mu)² ] / sum[w_i] )
 194  */
 195 void update_prior_sigma(float *mu, float *sigma, int num_players)
 196 {
 197         float nom = 0.0f, denom = 0.0f;
 198         for (int i = 0; i < num_players; ++i) {
 199                 for (unsigned j = 0; j < matches_for_player[i].size(); ++j) {
 200                         const match& m = matches_for_player[i][j];
 201
 202                         // Only count each match once.
 203                         if (m.other_player <= i) {
 204                                 continue;
 205                         }
 206
 207                         float mu1 = mu[i];
 208
 209                         float w = m.weight;
 210                         nom += w * ((mu1 - PRIOR_MU) * (mu1 - PRIOR_MU));
 211                         denom += w * 1.0f;
 212                 }
 213         }
 214
 215         prior_sigma = sqrt(nom / denom);
 216         if (!(prior_sigma > 40.0f)) {
 217                 prior_sigma = 40.0f;
 218         }
 219 }
 220
 221 float compute_logl(float z)
 222 {
 223         return -0.5 * (log(2.0f / M_PI) + z * z);
 224 }
 225
 226 float compute_total_logl(float *mu, float *sigma, int num_players)
 227 {
 228         float total_logl = 0.0f;
 229
 230         // Prior.
 231         for (int i = 0; i < num_players; ++i) {
 232                 total_logl += PRIOR_WEIGHT * compute_logl((mu[i] - PRIOR_MU) / prior_sigma);
 233         }
 234
 235         // Matches.
 236         for (int i = 0; i < num_players; ++i) {
 237                 for (unsigned j = 0; j < matches_for_player[i].size(); ++j) {
 238                         const match& m = matches_for_player[i][j];
 239
 240                         // Only count each match once.
 241                         if (m.other_player <= i) {
 242                                 continue;
 243                         }
 244
 245                         float mu1 = mu[i];
 246                         float mu2 = mu[m.other_player];
 247                         float sigma1 = sigma[i];
 248                         float sigma2 = sigma[m.other_player];
 249                         float sigma = sqrt(sigma1 * sigma1 + sigma2 * sigma2);
 250                         float mu = mu1 - mu2;
 251                         float x = m.margin;
 252                         float w = m.weight;
 253
 254                         total_logl += w * compute_logl((x - mu) / sigma);
 255                 }
 256         }
 257
 258         return total_logl;
 259 }
 260
 261 /*
 262  * Compute Hessian matrix of the negative log-likelihood, ie. for each term in logL:
 263  *
 264  * M_ij = D_i D_j (- logL) = -w / sigma²                                for i != j
 265  *                            w / sigma²                                for i == j
 266  *
 267  * Note that this does not depend on mu or the margin at all.
 268  */
 269 double hessian[MAX_PLAYERS][MAX_PLAYERS];
 270 void construct_hessian(const float *mu, const float *sigma, int num_players)
 271 {
 272         memset(hessian, 0, sizeof(hessian));
 273
 274         for (int i = 0; i < num_players; ++i) {
 275                 double sigma1 = sigma[i];
 276
 277                 for (unsigned k = 0; k < matches_for_player[i].size(); ++k) {
 278                         int j = matches_for_player[i][k].other_player;
 279
 280                         double sigma2 = sigma[j];
 281                         double sigma_sq = sigma1 * sigma1 + sigma2 * sigma2;
 282
 283                         float w = matches_for_player[i][k].weight;
 284
 285                         hessian[i][j] -= w / sigma_sq;
 286                         hessian[i][i] += w / sigma_sq;
 287                 }
 288         }
 289
 290         for (int i = 0; i < num_players; ++i) {
 291                 for (int j = 0; j < num_players; ++j) {
 292                         printf("%.12f ", hessian[i][j]);
 293                 }
 294                 printf("\n");
 295         }
 296 }
 297
 298 int main(int argc, char **argv)
 299 {
 300         int num_players;
 301         if (scanf("%d", &num_players) != 1) {
 302                 fprintf(stderr, "Could't read number of players\n");
 303                 exit(1);
 304         }
 305
 306         if (num_players > MAX_PLAYERS) {
 307                 fprintf(stderr, "Max %d players supported\n", MAX_PLAYERS);
 308                 exit(1);
 309         }
 310
 311         vector<string> players;
 312         map<string, int> player_map;
 313
 314         for (int i = 0; i < num_players; ++i) {
 315                 char buf[256];
 316                 if (scanf("%s", buf) != 1) {
 317                         fprintf(stderr, "Couldn't read player %d\n", i);
 318                         exit(1);
 319                 }
 320
 321                 players.push_back(buf);
 322                 player_map[buf] = i;
 323         }
 324
 325         int num_matches = 0;
 326         for ( ;; ) {
 327                 char pl1[256], pl2[256];
 328                 int score1, score2;
 329                 float weight;
 330
 331                 if (scanf("%s %s %d %d %f", pl1, pl2, &score1, &score2, &weight) != 5) {
 332                         //fprintf(stderr, "Read %d matches.\n", num_matches);
 333                         break;
 334                 }
 335
 336                 ++num_matches;
 337
 338                 if (player_map.count(pl1) == 0) {
 339                         fprintf(stderr, "Unknown player '%s'\n", pl1);
 340                         exit(1);
 341                 }
 342                 if (player_map.count(pl2) == 0) {
 343                         fprintf(stderr, "Unknown player '%s'\n", pl2);
 344                         exit(1);
 345                 }
 346
 347                 match m1;
 348                 m1.other_player = player_map[pl2];
 349                 m1.margin = score1 - score2;
 350                 m1.weight = weight;
 351                 matches_for_player[player_map[pl1]].push_back(m1);
 352
 353                 match m2;
 354                 m2.other_player = player_map[pl1];
 355                 m2.margin = score2 - score1;
 356                 m2.weight = weight;
 357                 matches_for_player[player_map[pl2]].push_back(m2);
 358         }
 359
 360         float mu[MAX_PLAYERS];
 361         float sigma[MAX_PLAYERS];
 362
 363         for (int i = 0; i < num_players; ++i) {
 364                 mu[i] = PRIOR_MU;
 365                 sigma[i] = 70.0f / sqrt(2.0f);
 366         }
 367
 368         for (int j = 0; j < 1000; ++j) {
 369                 float old_mu[MAX_PLAYERS];
 370                 float old_sigma[MAX_PLAYERS];
 371                 float old_prior_sigma = prior_sigma;
 372                 memcpy(old_mu, mu, sizeof(mu));
 373                 memcpy(old_sigma, sigma, sizeof(sigma));
 374                 for (int i = 0; i < num_players; ++i) {
 375                         update_mu(mu, sigma, i, matches_for_player[i]);
 376                 }
 377                 update_global_sigma(mu, sigma, num_players);
 378                 update_prior_sigma(mu, sigma, num_players);
 379                 /* for (int i = 0; i < num_players; ++i) {
 380                         update_sigma(mu, sigma, i, matches_for_player[i]);
 381                         dump_scores(players, mu, sigma, num_players);
 382                 } */
 383
 384                 float sumdiff = 0.0f;
 385                 for (int i = 0; i < num_players; ++i) {
 386                         sumdiff += (mu[i] - old_mu[i]) * (mu[i] - old_mu[i]);
 387                         sumdiff += (sigma[i] - old_sigma[i]) * (sigma[i] - old_sigma[i]);
 388                 }
 389                 sumdiff += (prior_sigma - old_prior_sigma) * (prior_sigma - old_prior_sigma);
 390                 if (sumdiff < EPSILON) {
 391                         //fprintf(stderr, "Converged after %d iterations. Stopping.\n", j);
 392                         printf("%d -1\n", j + 1);
 393                         break;
 394                 }
 395         }
 396
 397 #if DUMP_RAW
 398         dump_raw(mu, sigma, num_players);
 399 #else
 400         dump_scores(players, mu, sigma, num_players);
 401         //fprintf(stderr, "Optimal sigma: %f (two-player: %f)\n", sigma[0], sigma[0] * sqrt(2.0f));
 402         printf("%f -2\n", sigma[0]);
 403         printf("%f -3\n", prior_sigma);
 404
 405         float total_logl = compute_total_logl(mu, sigma, num_players);
 406         printf("%f -4\n", total_logl);
 407
 408 //      construct_hessian(mu, sigma, num_players);
 409 #endif
 410 }