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 };
  35 map<int, vector<match> > matches_for_player;
  36
  37 void dump_scores(const vector<string> &players, const float *mu, const float *sigma, int num_players)
  38 {
  39 #if 1
  40         for (int i = 0; i < num_players; ++i) {
  41                 fprintf(stderr, "%s=[%5.1f, %4.1f] ", players[i].c_str(), mu[i], sigma[i]);
  42         }
  43         fprintf(stderr, "\n");
  44 #else
  45         for (int i = 0; i < num_players; ++i) {
  46                 fprintf(stderr, "%5.1f ", mu[i]);
  47         }
  48         fprintf(stderr, "\n");
  49 #endif
  50 }
  51
  52 /*
  53  * diff(logL, mu1) = -(mu1 - mu2 - x) / sigma_c^2
  54  * maximizer for mu1 is given by: sum_i[ (1/sigma_c_i)^2 (mu1 - mu2_i - x_i) ] = 0
  55  *                                sum_i[ (1/sigma_c_i)^2 mu1 ] = sum_i [ (1/sigma_c_i)^2 ( mu2_i + x_i ) ]
  56  *                                mu1 = sum_i [ (1/sigma_c_i)^2 ( mu2_i + x_i ) ] / sum_i[ (1/sigma_c_i)^2 ]
  57  */
  58 void update_mu(float *mu, float *sigma, int player_num, const vector<match> &matches)
  59 {
  60         if (matches.empty()) {
  61                 return;
  62         }
  63
  64         float nom = 0.0f, denom = 0.0f;
  65         for (unsigned i = 0; i < matches.size(); ++i) {
  66                 float sigma1 = sigma[player_num];
  67                 float sigma2 = sigma[matches[i].other_player];
  68                 float inv_sigma_c2 = 1.0f / (sigma1 * sigma1 + sigma2 * sigma2);
  69                 float x = matches[i].margin; // / 70.0f;
  70
  71                 nom += (mu[matches[i].other_player] + x) * inv_sigma_c2;
  72                 denom += inv_sigma_c2;
  73         }
  74         mu[player_num] = nom / denom;
  75 }
  76
  77 /*
  78  * diff(logL, sigma1) = sigma1 (-sigma1² - sigma2² + (x - mu)²) / sigma_c²
  79  * maximizer for sigma1 is given by: sum_i[ (1/sigma_c_i)² sigma1 ((x - mu)² - (sigma1² + sigma2²) ] = 0
  80  *                                   sum_i[ (x - mu)² - sigma1² - sigma2² ] = 0                                  |: sigma1 != 0, sigma2 != 0
  81  *                                   sum_i[ (x - mu)² - sigma2² ] = sum[ sigma1² ]
  82  *                                   sigma1 = sqrt( sum_i[ (x - mu)² - sigma2² ] / N )
  83  */
  84 void update_sigma(float *mu, float *sigma, int player_num, const vector<match> &matches)
  85 {
  86         if (matches.size() < 2) {
  87                 return;
  88         }
  89
  90         float sum = 0.0f;
  91         for (unsigned i = 0; i < matches.size(); ++i) {
  92                 float mu1 = mu[player_num];
  93                 float mu2 = mu[matches[i].other_player];
  94                 float mu = mu1 - mu2;
  95                 float sigma2 = sigma[matches[i].other_player];
  96                 float x = matches[i].margin;
  97
  98                 //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);
  99                 sum += (x - mu) * (x - mu) - sigma2 * sigma2;
 100         }
 101
 102         if (sum <= 0) {
 103                 return;
 104         }
 105         //fprintf(stderr, "sum=%f\n", sum);
 106         sigma[player_num] = sqrt(sum / matches.size());
 107 }
 108
 109 void renormalize(float *mu, float *sigma, int num_players)
 110 {
 111         float avg = 0.0f;
 112         for (int i = 0; i < num_players; ++i) {
 113                 avg += mu[i];
 114         }
 115         float corr = 1500.0f - avg / num_players;
 116         for (int i = 0; i < num_players; ++i) {
 117                 mu[i] += corr;
 118         }
 119 }
 120
 121 /*
 122  * Compute Fisher information matrix of the log-likelihood, evaluated at the MLE,
 123 c
 124  * ie. M_ij = E[ (D_i logL) (D_j logL) ] = - sum( ( x - (mu_1 - mu_2) )² / sigma_c⁴ )  for i != j
 125  *                                       = - sum( 1 / sigma_c² )                     for i == j
 126  *
 127  * The Hessian matrix is generally zero and thus not as interesting.
 128  */
 129 void construct_fim(const float *mu, const float *sigma, int num_players)
 130 {
 131         float fim[MAX_PLAYERS][MAX_PLAYERS];
 132         memset(fim, 0, sizeof(fim));
 133
 134         for (int i = 0; i < num_players; ++i) {
 135                 float mu1 = mu[i];
 136                 float sigma1 = sigma[i];
 137
 138                 for (unsigned k = 0; k < matches_for_player[i].size(); ++k) {
 139                         int j = matches_for_player[i][k].other_player;
 140                         float mu2 = mu[j];
 141                         float sigma2 = sigma[j];
 142
 143                         float x = matches_for_player[i][k].margin;
 144                         float sigma_sq = sqrt(sigma1 * sigma1 + sigma2 * sigma2);
 145
 146                         fprintf(stderr, "exp_diff_sq=%f  sigma_sq=%f\n", (x - (mu1 - mu2)) * (x - (mu1 - mu2)), sigma_sq * sigma_sq);
 147
 148 #if 1
 149                         fim[i][i] += (x - (mu1 - mu2)) * (x - (mu1 - mu2)) / (sigma_sq * sigma_sq);
 150                         fim[i][j] -= (x - (mu1 - mu2)) * (x - (mu1 - mu2)) / (sigma_sq * sigma_sq);
 151 #else
 152                         fim[i][i] -= 1.0f / sigma_sq;
 153                         fim[i][j] += 1.0f / sigma_sq;
 154 #endif
 155                 }
 156
 157                 for (int j = 0; j < num_players; ++j) {
 158                         printf("%f ", fim[i][j]);
 159                 }
 160                 printf("\n");
 161         }
 162 }
 163
 164 int main(int argc, char **argv)
 165 {
 166         int num_players;
 167         if (scanf("%d", &num_players) != 1) {
 168                 fprintf(stderr, "Could't read number of players\n");
 169                 exit(1);
 170         }
 171
 172         if (num_players > MAX_PLAYERS) {
 173                 fprintf(stderr, "Max %d players supported\n", MAX_PLAYERS);
 174                 exit(1);
 175         }
 176
 177         vector<string> players;
 178         map<string, int> player_map;
 179
 180         for (int i = 0; i < num_players; ++i) {
 181                 char buf[256];
 182                 if (scanf("%s", buf) != 1) {
 183                         fprintf(stderr, "Couldn't read player %d\n", i);
 184                         exit(1);
 185                 }
 186
 187                 players.push_back(buf);
 188                 player_map[buf] = i;
 189         }
 190
 191         int num_matches = 0;
 192         for ( ;; ) {
 193                 char pl1[256], pl2[256];
 194                 int score1, score2;
 195
 196                 if (scanf("%s %s %d %d", pl1, pl2, &score1, &score2) != 4) {
 197                         fprintf(stderr, "Read %d matches.\n", num_matches);
 198                         break;
 199                 }
 200
 201                 ++num_matches;
 202
 203                 if (player_map.count(pl1) == 0) {
 204                         fprintf(stderr, "Unknown player '%s'\n", pl1);
 205                         exit(1);
 206                 }
 207                 if (player_map.count(pl2) == 0) {
 208                         fprintf(stderr, "Unknown player '%s'\n", pl2);
 209                         exit(1);
 210                 }
 211
 212                 match m1;
 213                 m1.other_player = player_map[pl2];
 214                 m1.margin = score1 - score2;
 215                 matches_for_player[player_map[pl1]].push_back(m1);
 216
 217                 match m2;
 218                 m2.other_player = player_map[pl1];
 219                 m2.margin = score2 - score1;
 220                 matches_for_player[player_map[pl2]].push_back(m2);
 221         }
 222
 223         float mu[MAX_PLAYERS];
 224         float sigma[MAX_PLAYERS];
 225
 226         for (int i = 0; i < num_players; ++i) {
 227                 mu[i] = 1500.0f;
 228                 sigma[i] = 70.0f / sqrt(2.0f);
 229         }
 230         renormalize(mu, sigma, num_players);
 231
 232         dump_scores(players, mu, sigma, num_players);
 233
 234         for (int j = 0; j < 100; ++j) {
 235                 float old_mu[MAX_PLAYERS];
 236                 float old_sigma[MAX_PLAYERS];
 237                 memcpy(old_mu, mu, sizeof(float) * MAX_PLAYERS);
 238                 memcpy(old_sigma, sigma, sizeof(float) * MAX_PLAYERS);
 239                 for (int i = 0; i < num_players; ++i) {
 240                         update_mu(mu, sigma, i, matches_for_player[i]);
 241                         renormalize(mu, sigma, num_players);
 242                         dump_scores(players, mu, sigma, num_players);
 243                 }
 244                 /* for (int i = 0; i < num_players; ++i) {
 245                         update_sigma(mu, sigma, i, matches_for_player[i]);
 246                         dump_scores(players, mu, sigma, num_players);
 247                 } */
 248                 bool any_difference = false;
 249                 for (int i = 0; i < num_players; ++i) {
 250                         if (fabs(mu[i] - old_mu[i]) > EPSILON ||
 251                             fabs(sigma[i] - old_sigma[i]) > EPSILON) {
 252                                 any_difference = true;
 253                                 break;
 254                         }
 255                 }
 256                 if (!any_difference) {
 257                         fprintf(stderr, "Converged after %d iterations. Stopping.\n", j);
 258                         break;
 259                 }
 260         }
 261
 262 //      construct_fim(mu, sigma, num_players);
 263 }