Invert the Hessian using Eigen instead of just using the diagonal. Gives slightly...

diff --git a/bayeswf.cpp b/bayeswf.cpp
index 71686c999cb944f2af6226d23b9fcf7e07d9115e..ca205ef47e4eb9cf75f9b27370d280e60d431b46 100644 (file)
--- a/bayeswf.cpp
+++ b/bayeswf.cpp
@@ -2,6 +2,8 @@
 #include <math.h>
 #include <string.h>
 #include <stdlib.h>
+#include <Eigen/Core>
+#include <Eigen/Eigenvalues>
 
 #include <map>
 #include <vector>
@@ -9,6 +11,7 @@
 #include <algorithm>
 
 using namespace std;
+using namespace Eigen;
 
 #define PRIOR_MU 1500
 #define PRIOR_WEIGHT 1.0
@@ -196,13 +199,13 @@ float compute_total_logl(float *mu, int num_players)
  *
  * Note that this does not depend on mu or the margin at all.
  */
-double hessian[MAX_PLAYERS][MAX_PLAYERS];
+double hessian[MAX_PLAYERS * MAX_PLAYERS];
 void construct_hessian(const float *mu, int num_players)
 {
        memset(hessian, 0, sizeof(hessian));
 
        for (int i = 0; i < num_players; ++i) {
-               hessian[i][i] += 1.0f / (prior_sigma * prior_sigma);
+               hessian[i * num_players + i] += 1.0f / (prior_sigma * prior_sigma);
        }
        for (unsigned i = 0; i < all_matches.size(); ++i) {
                const match &m = all_matches[i];
@@ -213,44 +216,29 @@ void construct_hessian(const float *mu, int num_players)
                double sigma_sq = global_sigma * global_sigma;
                float w = m.weight;
 
-               hessian[p1][p2] -= w / sigma_sq;
-               hessian[p2][p1] -= w / sigma_sq;
+               hessian[p1 * num_players + p2] -= w / sigma_sq;
+               hessian[p2 * num_players + p1] -= w / sigma_sq;
 
-               hessian[p1][p1] += w / sigma_sq;
-               hessian[p2][p2] += w / sigma_sq;
+               hessian[p1 * num_players + p1] += w / sigma_sq;
+               hessian[p2 * num_players + p2] += w / sigma_sq;
        }
 }
 
-// Compute uncertainty (stddev) of mu estimates, which is 1/sqrt(H_ii),
+// Compute uncertainty (stddev) of mu estimates, which is sqrt((H^-1)_ii),
 // where H is the Hessian (see construct_hessian()).
 void compute_mu_uncertainty(const float *mu, int num_players)
 {
-       memset(mu_stddev, 0, sizeof(mu_stddev));
-
-       // Temporarily use mu_stddev to store the diagonal of the Hessian.
-
-       // Prior.
+       Matrix<float, Dynamic, Dynamic> h(num_players, num_players);
        for (int i = 0; i < num_players; ++i) {
-               mu_stddev[i] += 1.0f / (prior_sigma * prior_sigma);
-       }
-
-       // Matches.
-       for (unsigned i = 0; i < all_matches.size(); ++i) {
-               const match &m = all_matches[i];
-
-               int p1 = m.player;
-               int p2 = m.other_player;
-
-               double sigma_sq = global_sigma * global_sigma;
-               float w = m.weight;
-
-               mu_stddev[p1] += w / sigma_sq;
-               mu_stddev[p2] += w / sigma_sq;
+               for (int j = 0; j < num_players; ++j) {
+                       h(i, j) = hessian[i * num_players + j];
+               }
        }
 
-       // Now convert to standard deviation.
+       // FIXME: Use pseudoinverse if applicable.
+       Matrix<float, Dynamic, Dynamic> ih = h.inverse();
        for (int i = 0; i < num_players; ++i) {
-               mu_stddev[i] = 1.0f / sqrt(mu_stddev[i]);
+               mu_stddev[i] = sqrt(ih(i, i));
        }
 }
 
@@ -357,6 +345,7 @@ int main(int argc, char **argv)
 #if DUMP_RAW
        dump_raw(mu, num_players);
 #else
+       construct_hessian(mu, num_players);
        compute_mu_uncertainty(mu, num_players);
        dump_scores(players, mu, mu_stddev, num_players);
        //fprintf(stderr, "Optimal sigma: %f (two-player: %f)\n", sigma[0], sigma[0] * sqrt(2.0f));
@@ -365,7 +354,5 @@ int main(int argc, char **argv)
 
        float total_logl = compute_total_logl(mu, num_players);
        printf("%f 0 -4\n", total_logl);
-
-//     construct_hessian(mu, sigma, num_players);
 #endif
 }

diff --git a/www/ratings-explained.html b/www/ratings-explained.html
index c2e5c7fd97d11504cdf3795325e6f0b140a4b64d..c45075b1ac8ab48a5338246e410c80da0c0e1cd1 100755 (executable)
--- a/www/ratings-explained.html
+++ b/www/ratings-explained.html
@@ -164,7 +164,7 @@
     <p>Dette var faktisk alt. Det skal sies at det sikkert er nok å ta tak i
       som ikke er blitt dekket her &ndash; for eksempel er det ikke beskrevet
       hvordan man regner ut <em>usikkerheten</em> i de estimerte ratingene
-      (hvilket er passe komplekst, og basert på en delvis invertering av
+      (hvilket er passe komplekst, og basert på å invertere
       <a href="http://en.wikipedia.org/wiki/Hessian_matrix">Hess-matrisen</a>
       til rimelighetsfunksjonen),
       eller hvordan modellen vekter kamper eldre kamper gis mindre betydning).</p>