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