X-Git-Url: https://git.sesse.net/?a=blobdiff_plain;f=bayeswf.cpp;h=cc591832a556f023ea6aa364a44ce5422a0b0740;hb=79ca58ab629b09e8fe59415578d2ee1ea5cb2079;hp=b7abb59dddf4e806588a3f9c12de9acc4bb1999e;hpb=156fbeb49e9aacf1b1d440b4ee352cd9bca457fc;p=wloh

diff --git a/bayeswf.cpp b/bayeswf.cpp
index b7abb59..cc59183 100644
--- 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,11 +11,17 @@
 #include <algorithm>
 
 using namespace std;
+using namespace Eigen;
 
+#define PRIOR_MU 500
+#define PRIOR_WEIGHT 1.0
 #define MAX_PLAYERS 4096
+#define DUMP_RAW 0
 
 float mu[MAX_PLAYERS];
-float sigma[MAX_PLAYERS];
+float mu_stddev[MAX_PLAYERS];
+float global_sigma = 70.0f;
+float prior_sigma = 70.0f;
 
 #define EPSILON 1e-3
 
@@ -29,13 +37,14 @@ float sigma[MAX_PLAYERS];
  */
 
 struct match {
-	int other_player;
+	int player, other_player;
 	int margin;
 	float weight;
 };
 map<int, vector<match> > matches_for_player;
+vector<match> all_matches;
 
-void dump_scores(const vector<string> &players, const float *mu, const float *sigma, int num_players)
+void dump_scores(const vector<string> &players, const float *mu, const float *mu_stddev, int num_players)
 {
 #if 0
 	for (int i = 0; i < num_players; ++i) {
@@ -49,9 +58,8 @@ void dump_scores(const vector<string> &players, const float *mu, const float *si
 	printf("\n");
 #else
 	for (int i = 0; i < num_players; ++i) {
-		printf("%5.1f %s\n", mu[i], players[i].c_str());
+		printf("%f %f %s\n", mu[i], mu_stddev[i], players[i].c_str());
 	}
-	printf("\n");
 #endif
 }
 
@@ -61,17 +69,24 @@ void dump_scores(const vector<string> &players, const float *mu, const float *si
  *                                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)
+void update_mu(float *mu, int player_num, const vector<match> &matches)
 {
 	if (matches.empty()) {
 		return;
 	}
 
 	float nom = 0.0f, denom = 0.0f;
+
+	// Prior.
+	{
+		float inv_sigma2 = 1.0f / (prior_sigma * prior_sigma);
+		nom += PRIOR_WEIGHT * PRIOR_MU * inv_sigma2;
+		denom += PRIOR_WEIGHT * inv_sigma2;
+	}
+
+	// All matches.
 	for (unsigned i = 0; i < matches.size(); ++i) {
-		float sigma1 = sigma[player_num];
-		float sigma2 = sigma[matches[i].other_player];
-		float inv_sigma_c2 = matches[i].weight / (sigma1 * sigma1 + sigma2 * sigma2);
+		float inv_sigma_c2 = matches[i].weight / (global_sigma * global_sigma);
 		float x = matches[i].margin; // / 70.0f;
 	
 		nom += (mu[matches[i].other_player] + x) * inv_sigma_c2;
@@ -80,84 +95,100 @@ void update_mu(float *mu, float *sigma, int player_num, const vector<match> &mat
 	mu[player_num] = nom / denom;
 }
 
+void dump_raw(const float *mu, int num_players)
+{
+	for (unsigned i = 0; i < all_matches.size(); ++i) {
+		const match& m = all_matches[i];
+
+		float mu1 = mu[m.player];
+		float mu2 = mu[m.other_player];
+		float sigma = global_sigma;
+		float mu = mu1 - mu2;
+		float x = m.margin;
+		float w = m.weight;
+
+		printf("%f %f\n", (x - mu) / sigma, w);
+	}
+}
+
 /*
- * diff(logL, sigma1) = sigma1 (-sigma1² - sigma2² + (x - mu)²) / sigma_c²
- * maximizer for sigma1 is given by: sum_i[ (1/sigma_c_i)² sigma1 ((x - mu)² - (sigma1² + sigma2²) ] = 0
- *                                   sum_i[ (x - mu)² - sigma1² - sigma2² ] = 0                                  |: sigma1 != 0, sigma2 != 0
- *                                   sum_i[ (x - mu)² - sigma2² ] = sum[ sigma1² ]
- *                                   sigma1 = sqrt( sum_i[ (x - mu)² - sigma2² ] / N )
+ * 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_sigma(float *mu, float *sigma, int player_num, const vector<match> &matches)
+void update_global_sigma(float *mu, int num_players)
 {
-	if (matches.size() < 2) {
-		return;
-	}
+	float nom = 0.0f, denom = 0.0f;
+	for (unsigned i = 0; i < all_matches.size(); ++i) {
+		const match& m = all_matches[i];
 
-	float sum = 0.0f;
-	for (unsigned i = 0; i < matches.size(); ++i) {
-		float mu1 = mu[player_num];
-		float mu2 = mu[matches[i].other_player];
+		float mu1 = mu[m.player];
+		float mu2 = mu[m.other_player];
 		float mu = mu1 - mu2;
-		float sigma2 = sigma[matches[i].other_player];
-		float x = matches[i].margin;
+		float x = m.margin;
+		float w = m.weight;
 
-		//fprintf(stderr, "x=%f mu=%f sigma2=%f   add %f-%f = %f\n", x, mu, sigma2, (x-mu)*(x-mu), sigma2*sigma2, (x - mu) * (x - mu) - sigma2 * sigma2);
-		sum += (x - mu) * (x - mu) - sigma2 * sigma2;
+		nom += w * ((x - mu) * (x - mu));
+		denom += w;
 	}
 
-	if (sum <= 0) {
-		return;
-	}
-	//fprintf(stderr, "sum=%f\n", sum);
-	sigma[player_num] = sqrt(sum / matches.size());
+	global_sigma = sqrt(nom / denom);
 }
 
 /*
- * diff(logL, sigma) = w ( (x - mu)² - sigma² ) / sigma³
+ * 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_global_sigma(float *mu, float *sigma, int num_players)
+void update_prior_sigma(float *mu, 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 mu1 = mu[i];
+
+		nom += ((mu1 - PRIOR_MU) * (mu1 - PRIOR_MU));
+		denom += 1.0f;
 	}
 
-	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;
+	prior_sigma = sqrt(nom / denom);
+	if (!(prior_sigma > 40.0f)) {
+		prior_sigma = 40.0f;
 	}
 }
 
-void renormalize(float *mu, float *sigma, int num_players)
+float compute_logl(float z)
+{
+	return -0.5 * (log(2.0f / M_PI) + z * z);
+}
+
+float compute_total_logl(float *mu, int num_players)
 {
-	float avg = 0.0f;
+	float total_logl = 0.0f;
+
+	// Prior.
 	for (int i = 0; i < num_players; ++i) {
-		avg += mu[i];
+		total_logl += PRIOR_WEIGHT * compute_logl((mu[i] - PRIOR_MU) / prior_sigma);
 	}
-	float corr = 1500.0f - avg / num_players;
-	for (int i = 0; i < num_players; ++i) {
-		mu[i] += corr;
+
+	// Matches.
+	for (unsigned i = 0; i < all_matches.size(); ++i) {
+		const match& m = all_matches[i];
+
+		float mu1 = mu[m.player];
+		float mu2 = mu[m.other_player];
+		float sigma = global_sigma;
+		float mu = mu1 - mu2;
+		float x = m.margin;
+		float w = m.weight;
+
+		total_logl += w * compute_logl((x - mu) / sigma);
 	}
+
+	return total_logl;
 }
 
 /*
@@ -168,32 +199,49 @@ void renormalize(float *mu, float *sigma, int num_players)
  *
  * Note that this does not depend on mu or the margin at all.
  */
-double hessian[MAX_PLAYERS][MAX_PLAYERS];
-void construct_hessian(const float *mu, const float *sigma, int num_players)
+Matrix<float, Dynamic, Dynamic> hessian;
+void construct_hessian(const float *mu, int num_players)
 {
-	memset(hessian, 0, sizeof(hessian));
+	hessian = Matrix<float, Dynamic, Dynamic>(num_players, num_players);
+	hessian.fill(0.0f);
 
 	for (int i = 0; i < num_players; ++i) {
-		double sigma1 = sigma[i];
+		hessian(i, i) += 1.0f / (prior_sigma * prior_sigma);
+	}
+	for (unsigned i = 0; i < all_matches.size(); ++i) {
+		const match &m = all_matches[i];
 
-		for (unsigned k = 0; k < matches_for_player[i].size(); ++k) {
-			int j = matches_for_player[i][k].other_player;
+		int p1 = m.player;
+		int p2 = m.other_player;
 
-			double sigma2 = sigma[j];
-			double sigma_sq = sigma1 * sigma1 + sigma2 * sigma2;
+		double sigma_sq = global_sigma * global_sigma;
+		float w = m.weight;
 
-			float w = matches_for_player[i][k].weight;
+		hessian(p1, p2) -= w / sigma_sq;
+		hessian(p2, p1) -= w / sigma_sq;
 
-			hessian[i][j] -= w / sigma_sq;
-			hessian[i][i] += w / sigma_sq;
-		}
+		hessian(p1, p1) += w / sigma_sq;
+		hessian(p2, p2) += w / sigma_sq;
 	}
+}
 
-	for (int i = 0; i < num_players; ++i) {
-		for (int j = 0; j < num_players; ++j) {
-			printf("%.12f ", hessian[i][j]);
+// 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, const vector<string> &players)
+{
+	// FIXME: Use pseudoinverse if applicable.
+	Matrix<float, Dynamic, Dynamic> ih = hessian.inverse();
+	for (unsigned i = 0; i < players.size(); ++i) {
+		mu_stddev[i] = sqrt(ih(i, i));
+	}
+
+	for (unsigned i = 0; i < players.size(); ++i) {
+		for (unsigned j = 0; j < players.size(); ++j) {
+			printf("covariance %s %s %f\n",
+			       players[i].c_str(),
+			       players[j].c_str(),
+			       ih(i, j));
 		}
-		printf("\n");
 	}
 }
 
@@ -231,7 +279,7 @@ int main(int argc, char **argv)
 		float weight;
 
 		if (scanf("%s %s %d %d %f", pl1, pl2, &score1, &score2, &weight) != 5) {
-			fprintf(stderr, "Read %d matches.\n", num_matches);
+			//fprintf(stderr, "Read %d matches.\n", num_matches);
 			break;
 		}
 
@@ -247,37 +295,38 @@ int main(int argc, char **argv)
 		}
 
 		match m1;
+		m1.player = player_map[pl1];
 		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.player = player_map[pl2];
 		m2.other_player = player_map[pl1];
 		m2.margin = score2 - score1;
 		m2.weight = weight;
 		matches_for_player[player_map[pl2]].push_back(m2);
+
+		all_matches.push_back(m1);
 	}
 
 	float mu[MAX_PLAYERS];
-	float sigma[MAX_PLAYERS];
 
 	for (int i = 0; i < num_players; ++i) {
-		mu[i] = 1500.0f;
-		sigma[i] = 70.0f / sqrt(2.0f);
+		mu[i] = PRIOR_MU;
 	}
-	renormalize(mu, sigma, num_players);
 
 	for (int j = 0; j < 1000; ++j) {
 		float old_mu[MAX_PLAYERS];
-		float old_sigma[MAX_PLAYERS];
+		float old_global_sigma = global_sigma;
+		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_mu(mu, i, matches_for_player[i]);
 		}
-		update_global_sigma(mu, sigma, num_players);
+		update_global_sigma(mu, num_players);
+		update_prior_sigma(mu, 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);
@@ -286,15 +335,27 @@ int main(int argc, char **argv)
 		float sumdiff = 0.0f;
 		for (int i = 0; i < num_players; ++i) {
 			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);
+		sumdiff += (global_sigma - old_global_sigma) * (global_sigma - old_global_sigma);
 		if (sumdiff < EPSILON) {
-			fprintf(stderr, "Converged after %d iterations. Stopping.\n", j);
+			//fprintf(stderr, "Converged after %d iterations. Stopping.\n", j);
+			printf("aux_param num_iterations %d\n", j + 1);
 			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_hessian(mu, sigma, num_players);
+#if DUMP_RAW
+	dump_raw(mu, num_players);
+#else
+	construct_hessian(mu, num_players);
+	compute_mu_uncertainty(mu, 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));
+	printf("aux_param score_stddev %f\n", global_sigma / sqrt(2.0f));
+	printf("aux_param rating_prior_stddev %f\n", prior_sigma);
+
+	float total_logl = compute_total_logl(mu, num_players);
+	printf("aux_param total_log_likelihood %f\n", total_logl);
+#endif
 }