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