Compute the Hessian instead of the Fisher information matrix. The numerics here suck.

diff --git a/bayeswf.cpp b/bayeswf.cpp
index 48d52325620c0542b4f61553e2d0d075d906371c..56214b5494b5647fb718c63040ef29652fb5d4ca 100644 (file)
--- a/bayeswf.cpp
+++ b/bayeswf.cpp
@@ -161,43 +161,35 @@ 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:
  *
- * The Hessian matrix is generally zero and thus not as interesting.
+ * M_ij = D_i D_j (- logL) = -1 / sigma²                                for i != j
+ *                            1 / sigma²                                for i == j
+ *
+ * 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);
 
-                       fprintf(stderr, "exp_diff_sq=%f  sigma_sq=%f\n", (x - (mu1 - mu2)) * (x - (mu1 - mu2)), sigma_sq * sigma_sq);
+                       double sigma2 = sigma[j];
+                       double sigma_sq = sigma1 * sigma1 + sigma2 * sigma2;
 
-#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] -= 1.0f / sigma_sq;
+                       hessian[i][i] += 1.0f / 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");
        }
@@ -302,5 +294,5 @@ int main(int argc, char **argv)
        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);
 }