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