X-Git-Url: https://git.sesse.net/?a=blobdiff_plain;f=bayeswf.cpp;h=cc591832a556f023ea6aa364a44ce5422a0b0740;hb=1fac0677bbcba5dec26e11aa661c97e38f6c3d40;hp=9365a1068ea460efb41e0008ef8543e0e5dd9053;hpb=85c15ca0ba6e416371abb36b0c1c63a4a02e64fd;p=wloh diff --git a/bayeswf.cpp b/bayeswf.cpp index 9365a10..cc59183 100644 --- a/bayeswf.cpp +++ b/bayeswf.cpp @@ -2,6 +2,8 @@ #include #include #include +#include +#include #include #include @@ -9,14 +11,16 @@ #include 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 @@ -40,7 +44,7 @@ struct match { map > matches_for_player; vector all_matches; -void dump_scores(const vector &players, const float *mu, const float *sigma, int num_players) +void dump_scores(const vector &players, const float *mu, const float *mu_stddev, int num_players) { #if 0 for (int i = 0; i < num_players; ++i) { @@ -54,7 +58,7 @@ void dump_scores(const vector &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 } @@ -65,7 +69,7 @@ void dump_scores(const vector &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 &matches) +void update_mu(float *mu, int player_num, const vector &matches) { if (matches.empty()) { return; @@ -82,9 +86,7 @@ void update_mu(float *mu, float *sigma, int player_num, const vector &mat // 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; @@ -93,16 +95,14 @@ void update_mu(float *mu, float *sigma, int player_num, const vector &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 (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 sigma1 = sigma[m.player]; - float sigma2 = sigma[m.other_player]; - float sigma = sqrt(sigma1 * sigma1 + sigma2 * sigma2); + float sigma = global_sigma; float mu = mu1 - mu2; float x = m.margin; float w = m.weight; @@ -111,38 +111,6 @@ void dump_raw(const float *mu, const float *sigma, int num_players) } } -/* - * 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 &matches) -{ - if (matches.size() < 2) { - return; - } - - 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 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; - } - - if (sum <= 0) { - return; - } - //fprintf(stderr, "sum=%f\n", sum); - sigma[player_num] = sqrt(sum / matches.size()); -} - /* * diff(logL, sigma) = w ( (x - mu)² - sigma² ) / sigma³ * maximizer for sigma is given by: sum_i[ (w_i/sigma)³ ((x - mu)² - sigma²) ] = 0 @@ -150,7 +118,7 @@ void update_sigma(float *mu, float *sigma, int player_num, const vector & * 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 (unsigned i = 0; i < all_matches.size(); ++i) { @@ -166,10 +134,7 @@ void update_global_sigma(float *mu, float *sigma, int num_players) denom += w; } - 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; - } + global_sigma = sqrt(nom / denom); } /* @@ -179,7 +144,7 @@ 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) { @@ -200,7 +165,7 @@ float compute_logl(float z) return -0.5 * (log(2.0f / M_PI) + z * z); } -float compute_total_logl(float *mu, float *sigma, int num_players) +float compute_total_logl(float *mu, int num_players) { float total_logl = 0.0f; @@ -215,9 +180,7 @@ float compute_total_logl(float *mu, float *sigma, int num_players) float mu1 = mu[m.player]; float mu2 = mu[m.other_player]; - float sigma1 = sigma[m.player]; - float sigma2 = sigma[m.other_player]; - float sigma = sqrt(sigma1 * sigma1 + sigma2 * sigma2); + float sigma = global_sigma; float mu = mu1 - mu2; float x = m.margin; float w = m.weight; @@ -236,32 +199,49 @@ float compute_total_logl(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 hessian; +void construct_hessian(const float *mu, int num_players) { - memset(hessian, 0, sizeof(hessian)); + hessian = Matrix(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 &players) +{ + // FIXME: Use pseudoinverse if applicable. + Matrix 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"); } } @@ -332,24 +312,21 @@ int main(int argc, char **argv) } float mu[MAX_PLAYERS]; - float sigma[MAX_PLAYERS]; for (int i = 0; i < num_players; ++i) { mu[i] = PRIOR_MU; - sigma[i] = 70.0f / sqrt(2.0f); } 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); @@ -358,27 +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); - - float total_logl = compute_total_logl(mu, sigma, num_players); - printf("%f -4\n", total_logl); + 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 }