]> git.sesse.net Git - wloh/blobdiff - bayeswf.cpp
Take weight into account in Hessian.
[wloh] / bayeswf.cpp
index 8c6f3b061f23da7bf74e0f1768d7a070d54ba01a..b7abb59dddf4e806588a3f9c12de9acc4bb1999e 100644 (file)
@@ -31,6 +31,7 @@ float sigma[MAX_PLAYERS];
 struct match {
        int other_player;
        int margin;
+       float weight;
 };
 map<int, vector<match> > matches_for_player;
 
@@ -55,10 +56,10 @@ void dump_scores(const vector<string> &players, const float *mu, const float *si
 }
 
 /*
- * diff(logL, mu1) = -(mu1 - mu2 - x) / sigma_c^2
- * maximizer for mu1 is given by: sum_i[ (1/sigma_c_i)^2 (mu1 - mu2_i - x_i) ] = 0
- *                                sum_i[ (1/sigma_c_i)^2 mu1 ] = sum_i [ (1/sigma_c_i)^2 ( mu2_i + x_i ) ]
- *                                mu1 = sum_i [ (1/sigma_c_i)^2 ( mu2_i + x_i ) ] / sum_i[ (1/sigma_c_i)^2 ]
+ * diff(logL, mu1) = -w * (mu1 - mu2 - x) / sigma_c^2
+ * maximizer for mu1 is given by: sum_i[ (w_i/sigma_c_i)^2 (mu1 - mu2_i - x_i) ] = 0
+ *                                sum_i[ (w_i/sigma_c_i)^2 mu1 ] = sum_i [ (w_i/sigma_c_i)^2 ( mu2_i + x_i ) ]
+ *                                mu1 = sum_i [ (w_i/sigma_c_i)^2 ( mu2_i + x_i ) ] / sum_i[ (w_i/sigma_c_i)^2 ]
  */
 void update_mu(float *mu, float *sigma, int player_num, const vector<match> &matches)
 {
@@ -70,7 +71,7 @@ void update_mu(float *mu, float *sigma, int player_num, const vector<match> &mat
        for (unsigned i = 0; i < matches.size(); ++i) {
                float sigma1 = sigma[player_num];
                float sigma2 = sigma[matches[i].other_player];
-               float inv_sigma_c2 = 1.0f / (sigma1 * sigma1 + sigma2 * sigma2);
+               float inv_sigma_c2 = matches[i].weight / (sigma1 * sigma1 + sigma2 * sigma2);
                float x = matches[i].margin; // / 70.0f;
        
                nom += (mu[matches[i].other_player] + x) * inv_sigma_c2;
@@ -111,6 +112,42 @@ void update_sigma(float *mu, float *sigma, int player_num, const vector<match> &
        sigma[player_num] = sqrt(sum / matches.size());
 }
 
+/*
+ * diff(logL, sigma) = w ( (x - mu)² - sigma² ) / sigma³
+ * maximizer for sigma is given by: sum_i[ (w_i/sigma)³ ((x - mu)² - sigma²) ] = 0
+ *                                   sum_i[ w_i ( (x - mu)² - sigma² ) ] = 0                            |: sigma != 0
+ *                                   sum_i[ w_i (x - mu)² ] = sum[ w_i sigma² ]
+ *                                   sigma = sqrt( sum_i[ w_i (x - mu)² ] / sum[w_i] )
+ */
+void update_global_sigma(float *mu, float *sigma, int num_players)
+{
+       float nom = 0.0f, denom = 0.0f;
+       for (int i = 0; i < num_players; ++i) {
+               for (unsigned j = 0; j < matches_for_player[i].size(); ++j) {
+                       const match& m = matches_for_player[i][j];
+
+                       // Only count each match once.
+                       if (m.other_player <= i) {
+                               continue;
+                       }
+
+                       float mu1 = mu[i];
+                       float mu2 = mu[m.other_player];
+                       float mu = mu1 - mu2;
+                       float x = m.margin;
+                       float w = m.weight;
+
+                       nom += w * ((x - mu) * (x - mu));
+                       denom += w;
+               }
+       }
+
+       float best_sigma = sqrt(nom / denom) / sqrt(2.0f);  // Divide evenly between the two players.
+       for (int i = 0; i < num_players; ++i) {
+               sigma[i] = best_sigma;
+       }
+}
+
 void renormalize(float *mu, float *sigma, int num_players)
 {
        float avg = 0.0f;
@@ -124,43 +161,37 @@ void renormalize(float *mu, float *sigma, int num_players)
 }
 
 /*
- * Compute Fisher information matrix of the log-likelihood, evaluated at the MLE,
-c
- * ie. M_ij = E[ (D_i logL) (D_j logL) ] = - sum( ( x - (mu_1 - mu_2) )² / sigma_c⁴ )  for i != j
- *                                       = - sum( 1 / sigma_c² )                     for i == j
+ * Compute Hessian matrix of the negative log-likelihood, ie. for each term in logL:
+ *
+ * M_ij = D_i D_j (- logL) = -w / sigma²                                for i != j
+ *                            w / sigma²                                for i == j
  *
- * The Hessian matrix is generally zero and thus not as interesting.
+ * Note that this does not depend on mu or the margin at all.
  */
-void construct_fim(const float *mu, const float *sigma, int num_players)
+double hessian[MAX_PLAYERS][MAX_PLAYERS];
+void construct_hessian(const float *mu, const float *sigma, int num_players)
 {
-       float fim[MAX_PLAYERS][MAX_PLAYERS];
-       memset(fim, 0, sizeof(fim));
+       memset(hessian, 0, sizeof(hessian));
 
        for (int i = 0; i < num_players; ++i) {
-               float mu1 = mu[i];
-               float sigma1 = sigma[i];
+               double sigma1 = sigma[i];
 
                for (unsigned k = 0; k < matches_for_player[i].size(); ++k) {
                        int j = matches_for_player[i][k].other_player;
-                       float mu2 = mu[j];
-                       float sigma2 = sigma[j];
 
-                       float x = matches_for_player[i][k].margin;
-                       float sigma_sq = sqrt(sigma1 * sigma1 + sigma2 * sigma2);
+                       double sigma2 = sigma[j];
+                       double sigma_sq = sigma1 * sigma1 + sigma2 * sigma2;
 
-                       fprintf(stderr, "exp_diff_sq=%f  sigma_sq=%f\n", (x - (mu1 - mu2)) * (x - (mu1 - mu2)), sigma_sq * sigma_sq);
+                       float w = matches_for_player[i][k].weight;
 
-#if 1
-                       fim[i][i] += (x - (mu1 - mu2)) * (x - (mu1 - mu2)) / (sigma_sq * sigma_sq);
-                       fim[i][j] -= (x - (mu1 - mu2)) * (x - (mu1 - mu2)) / (sigma_sq * sigma_sq);
-#else
-                       fim[i][i] -= 1.0f / sigma_sq;
-                       fim[i][j] += 1.0f / sigma_sq;
-#endif
+                       hessian[i][j] -= w / sigma_sq;
+                       hessian[i][i] += w / sigma_sq;
                }
+       }
 
+       for (int i = 0; i < num_players; ++i) {
                for (int j = 0; j < num_players; ++j) {
-                       printf("%f ", fim[i][j]);
+                       printf("%.12f ", hessian[i][j]);
                }
                printf("\n");
        }
@@ -197,8 +228,9 @@ int main(int argc, char **argv)
        for ( ;; ) {
                char pl1[256], pl2[256];
                int score1, score2;
+               float weight;
 
-               if (scanf("%s %s %d %d", pl1, pl2, &score1, &score2) != 4) {
+               if (scanf("%s %s %d %d %f", pl1, pl2, &score1, &score2, &weight) != 5) {
                        fprintf(stderr, "Read %d matches.\n", num_matches);
                        break;
                }
@@ -217,11 +249,13 @@ int main(int argc, char **argv)
                match m1;
                m1.other_player = player_map[pl2];
                m1.margin = score1 - score2;
+               m1.weight = weight;
                matches_for_player[player_map[pl1]].push_back(m1);
 
                match m2;
                m2.other_player = player_map[pl1];
                m2.margin = score2 - score1;
+               m2.weight = weight;
                matches_for_player[player_map[pl2]].push_back(m2);
        }
 
@@ -234,33 +268,33 @@ int main(int argc, char **argv)
        }
        renormalize(mu, sigma, num_players);
 
-       for (int j = 0; j < 100; ++j) {
+       for (int j = 0; j < 1000; ++j) {
                float old_mu[MAX_PLAYERS];
                float old_sigma[MAX_PLAYERS];
-               memcpy(old_mu, mu, sizeof(float) * MAX_PLAYERS);
-               memcpy(old_sigma, sigma, sizeof(float) * MAX_PLAYERS);
+               memcpy(old_mu, mu, sizeof(mu));
+               memcpy(old_sigma, sigma, sizeof(sigma));
                for (int i = 0; i < num_players; ++i) {
                        update_mu(mu, sigma, i, matches_for_player[i]);
                        renormalize(mu, sigma, num_players);
                }
+               update_global_sigma(mu, sigma, num_players);
                /* for (int i = 0; i < num_players; ++i) {
                        update_sigma(mu, sigma, i, matches_for_player[i]);
                        dump_scores(players, mu, sigma, num_players);
                } */
-               bool any_difference = false;
+
+               float sumdiff = 0.0f;
                for (int i = 0; i < num_players; ++i) {
-                       if (fabs(mu[i] - old_mu[i]) > EPSILON ||
-                           fabs(sigma[i] - old_sigma[i]) > EPSILON) {
-                               any_difference = true;
-                               break;
-                       }
+                       sumdiff += (mu[i] - old_mu[i]) * (mu[i] - old_mu[i]);
+                       sumdiff += (sigma[i] - old_sigma[i]) * (sigma[i] - old_sigma[i]);
                }
-               if (!any_difference) {
+               if (sumdiff < EPSILON) {
                        fprintf(stderr, "Converged after %d iterations. Stopping.\n", j);
                        break;
                }
        }
        dump_scores(players, mu, sigma, num_players);
+       fprintf(stderr, "Optimal sigma: %f (two-player: %f)\n", sigma[0], sigma[0] * sqrt(2.0f));
 
-//     construct_fim(mu, sigma, num_players);
+       construct_hessian(mu, sigma, num_players);
 }