X-Git-Url: https://git.sesse.net/?a=blobdiff_plain;f=bayeswf.cpp;h=cc591832a556f023ea6aa364a44ce5422a0b0740;hb=1fac0677bbcba5dec26e11aa661c97e38f6c3d40;hp=c54f72ab03d82cb7edc5c18f209e7620a3490ca1;hpb=7a21b8634e3bb28e817746782db718597ac0d7ea;p=wloh

diff --git a/bayeswf.cpp b/bayeswf.cpp
index c54f72a..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,13 +11,16 @@
 #include <algorithm>
 
 using namespace std;
+using namespace Eigen;
 
-#define PRIOR_MU 1500
+#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
@@ -32,13 +37,14 @@ float prior_sigma = 70.0f;
  */
 
 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) {
@@ -52,7 +58,7 @@ 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("%f %s\n", mu[i], players[i].c_str());
+		printf("%f %f %s\n", mu[i], mu_stddev[i], players[i].c_str());
 	}
 #endif
 }
@@ -63,7 +69,7 @@ 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;
@@ -74,15 +80,13 @@ void update_mu(float *mu, float *sigma, int player_num, const vector<match> &mat
 	// Prior.
 	{
 		float inv_sigma2 = 1.0f / (prior_sigma * prior_sigma);
-		nom += PRIOR_MU * inv_sigma2;
-		denom += inv_sigma2;
+		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;
@@ -91,61 +95,20 @@ 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, const float *sigma, int num_players)
+void dump_raw(const float *mu, int num_players)
 {
-	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 sigma1 = sigma[i];
-			float sigma2 = sigma[m.other_player];
-			float sigma = sqrt(sigma1 * sigma1 + sigma2 * sigma2);
-			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 )
- */
-void update_sigma(float *mu, float *sigma, int player_num, const vector<match> &matches)
-{
-	if (matches.size() < 2) {
-		return;
-	}
+	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 sigma = global_sigma;
 		float mu = mu1 - mu2;
-		float sigma2 = sigma[matches[i].other_player];
-		float x = matches[i].margin;
-
-		//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;
-	}
+		float x = m.margin;
+		float w = m.weight;
 
-	if (sum <= 0) {
-		return;
+		printf("%f %f\n", (x - mu) / sigma, w);
 	}
-	//fprintf(stderr, "sum=%f\n", sum);
-	sigma[player_num] = sqrt(sum / matches.size());
 }
 
 /*
@@ -155,33 +118,23 @@ void update_sigma(float *mu, float *sigma, int player_num, const vector<match> &
  *                                   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_global_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;
-		}
-	}
+	for (unsigned i = 0; i < all_matches.size(); ++i) {
+		const match& m = all_matches[i];
 
-	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;
+		float mu1 = mu[m.player];
+		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;
 	}
+
+	global_sigma = sqrt(nom / denom);
 }
 
 /*
@@ -191,27 +144,51 @@ void update_global_sigma(float *mu, float *sigma, int num_players)
  *                                   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)
+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];
+		float mu1 = mu[i];
 
-			// Only count each match once.
-			if (m.other_player <= i) {
-				continue;
-			}
+		nom += ((mu1 - PRIOR_MU) * (mu1 - PRIOR_MU));
+		denom += 1.0f;
+	}
 
-			float mu1 = mu[i];
+	prior_sigma = sqrt(nom / denom);
+	if (!(prior_sigma > 40.0f)) {
+		prior_sigma = 40.0f;
+	}
+}
 
-			float w = m.weight;
-			nom += w * ((mu1 - PRIOR_MU) * (mu1 - PRIOR_MU));
-			denom += w * 1.0f;
-		}
+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 total_logl = 0.0f;
+
+	// Prior.
+	for (int i = 0; i < num_players; ++i) {
+		total_logl += PRIOR_WEIGHT * compute_logl((mu[i] - PRIOR_MU) / prior_sigma);
 	}
 
-	prior_sigma = sqrt(nom / denom);
+	// 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;
 }
 
 /*
@@ -222,32 +199,49 @@ void update_prior_sigma(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");
 	}
 }
 
@@ -301,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;
 	}
 
 	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]);
+			update_mu(mu, i, matches_for_player[i]);
 		}
-		update_global_sigma(mu, sigma, num_players);
-		update_prior_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);
@@ -340,24 +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);
-			printf("%d -1\n", j + 1);
+			printf("aux_param num_iterations %d\n", j + 1);
 			break;
 		}
 	}
 
 #if DUMP_RAW
-	dump_raw(mu, sigma, num_players);
+	dump_raw(mu, num_players);
 #else
-	dump_scores(players, mu, sigma, num_players);
+	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("%f -2\n", sigma[0]);
-	printf("%f -3\n", prior_sigma);
+	printf("aux_param score_stddev %f\n", global_sigma / sqrt(2.0f));
+	printf("aux_param rating_prior_stddev %f\n", prior_sigma);
 
-//	construct_hessian(mu, sigma, num_players);
+	float total_logl = compute_total_logl(mu, num_players);
+	printf("aux_param total_log_likelihood %f\n", total_logl);
 #endif
 }