From b0f2024b2a7e87cbe8d59f4c0365e0b49e16f532 Mon Sep 17 00:00:00 2001
From: "Steinar H. Gunderson" <Steinar H. Gunderson sesse@debian.org>
Date: Sat, 17 Mar 2012 14:10:40 +0100
Subject: [PATCH] Optimize for the prior sigma, and stop renormalizing since
 the prior will take care of that.

---
 bayeswf.cpp | 41 ++++++++++++++++++++++++++++++-----------
 1 file changed, 30 insertions(+), 11 deletions(-)

diff --git a/bayeswf.cpp b/bayeswf.cpp
index 176d345..5eba479 100644
--- a/bayeswf.cpp
+++ b/bayeswf.cpp
@@ -11,11 +11,11 @@
 using namespace std;
 
 #define PRIOR_MU 1500
-#define PRIOR_SIGMA 100
 #define MAX_PLAYERS 4096
 
 float mu[MAX_PLAYERS];
 float sigma[MAX_PLAYERS];
+float prior_sigma = 70.0f;
 
 #define EPSILON 1e-3
 
@@ -72,7 +72,7 @@ void update_mu(float *mu, float *sigma, int player_num, const vector<match> &mat
 
 	// Prior.
 	{
-		float inv_sigma2 = 1.0f / (PRIOR_SIGMA * PRIOR_SIGMA);
+		float inv_sigma2 = 1.0f / (prior_sigma * prior_sigma);
 		nom += PRIOR_MU * inv_sigma2;
 		denom += inv_sigma2;
 	}
@@ -158,16 +158,33 @@ void update_global_sigma(float *mu, float *sigma, int num_players)
 	}
 }
 
-void renormalize(float *mu, float *sigma, int num_players)
+/*
+ * diff(priorlogL, 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_prior_sigma(float *mu, float *sigma, int num_players)
 {
-	float avg = 0.0f;
-	for (int i = 0; i < num_players; ++i) {
-		avg += mu[i];
-	}
-	float corr = 1500.0f - avg / num_players;
+	float nom = 0.0f, denom = 0.0f;
 	for (int i = 0; i < num_players; ++i) {
-		mu[i] += corr;
+		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];
+
+			nom += ((mu1 - PRIOR_MU) * (mu1 - PRIOR_MU));
+			denom += 1.0f;
+		}
 	}
+
+	prior_sigma = sqrt(nom / denom);
 }
 
 /*
@@ -276,18 +293,18 @@ int main(int argc, char **argv)
 		mu[i] = 1500.0f;
 		sigma[i] = 70.0f / sqrt(2.0f);
 	}
-	renormalize(mu, sigma, num_players);
 
 	for (int j = 0; j < 1000; ++j) {
 		float old_mu[MAX_PLAYERS];
 		float old_sigma[MAX_PLAYERS];
+		float old_prior_sigma = prior_sigma;
 		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);
+		update_prior_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);
@@ -298,6 +315,7 @@ int main(int argc, char **argv)
 			sumdiff += (mu[i] - old_mu[i]) * (mu[i] - old_mu[i]);
 			sumdiff += (sigma[i] - old_sigma[i]) * (sigma[i] - old_sigma[i]);
 		}
+		sumdiff += (prior_sigma - old_prior_sigma) * (prior_sigma - old_prior_sigma);
 		if (sumdiff < EPSILON) {
 			//fprintf(stderr, "Converged after %d iterations. Stopping.\n", j);
 			printf("%d -1\n", j + 1);
@@ -307,6 +325,7 @@ 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));
 	printf("%f -2\n", sigma[0]);
+	printf("%f -3\n", prior_sigma);
 
 //	construct_hessian(mu, sigma, num_players);
 }
-- 
2.39.5