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