From 261a1578a49b9fa228d40b48af04dd26aed6b60f Mon Sep 17 00:00:00 2001 From: "Steinar H. Gunderson" Date: Sat, 17 Mar 2012 02:29:41 +0100 Subject: [PATCH] Compute the Hessian instead of the Fisher information matrix. The numerics here suck. --- bayeswf.cpp | 42 +++++++++++++++++------------------------- 1 file changed, 17 insertions(+), 25 deletions(-) diff --git a/bayeswf.cpp b/bayeswf.cpp index 48d5232..56214b5 100644 --- 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); } -- 2.39.2