X-Git-Url: https://git.sesse.net/?a=blobdiff_plain;f=bayeswf.cpp;h=cc591832a556f023ea6aa364a44ce5422a0b0740;hb=79ca58ab629b09e8fe59415578d2ee1ea5cb2079;hp=8c6f3b061f23da7bf74e0f1768d7a070d54ba01a;hpb=aac670b21f109945d5766dcb15946a2ea44da807;p=wloh diff --git a/bayeswf.cpp b/bayeswf.cpp index 8c6f3b0..cc59183 100644 --- a/bayeswf.cpp +++ b/bayeswf.cpp @@ -2,6 +2,8 @@ #include #include #include +#include +#include #include #include @@ -9,11 +11,17 @@ #include 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,12 +37,14 @@ float sigma[MAX_PLAYERS]; */ struct match { - int other_player; + int player, other_player; int margin; + float weight; }; 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) { @@ -48,29 +58,35 @@ void dump_scores(const vector &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 } /* - * diff(logL, mu1) = -(mu1 - mu2 - x) / sigma_c^2 - * maximizer for mu1 is given by: sum_i[ (1/sigma_c_i)^2 (mu1 - mu2_i - x_i) ] = 0 - * sum_i[ (1/sigma_c_i)^2 mu1 ] = sum_i [ (1/sigma_c_i)^2 ( mu2_i + x_i ) ] - * mu1 = sum_i [ (1/sigma_c_i)^2 ( mu2_i + x_i ) ] / sum_i[ (1/sigma_c_i)^2 ] + * diff(logL, mu1) = -w * (mu1 - mu2 - x) / sigma_c^2 + * maximizer for mu1 is given by: sum_i[ (w_i/sigma_c_i)^2 (mu1 - mu2_i - x_i) ] = 0 + * 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; } 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 = 1.0f / (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; @@ -79,90 +95,153 @@ void update_mu(float *mu, float *sigma, int player_num, const vector &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 &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); } -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, int num_players) { - float avg = 0.0f; + float nom = 0.0f, denom = 0.0f; for (int i = 0; i < num_players; ++i) { - avg += mu[i]; + float mu1 = mu[i]; + + nom += ((mu1 - PRIOR_MU) * (mu1 - PRIOR_MU)); + denom += 1.0f; + } + + prior_sigma = sqrt(nom / denom); + if (!(prior_sigma > 40.0f)) { + prior_sigma = 40.0f; } - float corr = 1500.0f - avg / 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 total_logl = 0.0f; + + // Prior. for (int i = 0; i < num_players; ++i) { - mu[i] += corr; + total_logl += PRIOR_WEIGHT * compute_logl((mu[i] - PRIOR_MU) / prior_sigma); } + + // 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; } /* - * 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: + * + * M_ij = D_i D_j (- logL) = -w / sigma² for i != j + * w / sigma² for i == j * - * The Hessian matrix is generally zero and thus not as interesting. + * 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) +Matrix hessian; +void construct_hessian(const float *mu, int num_players) { - float fim[MAX_PLAYERS][MAX_PLAYERS]; - memset(fim, 0, sizeof(fim)); + hessian = Matrix(num_players, num_players); + hessian.fill(0.0f); for (int i = 0; i < num_players; ++i) { - float mu1 = mu[i]; - float 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; - float mu2 = mu[j]; - float sigma2 = sigma[j]; + int p1 = m.player; + int p2 = m.other_player; - float x = matches_for_player[i][k].margin; - float sigma_sq = sqrt(sigma1 * sigma1 + sigma2 * sigma2); + double sigma_sq = global_sigma * global_sigma; + float w = m.weight; - fprintf(stderr, "exp_diff_sq=%f sigma_sq=%f\n", (x - (mu1 - mu2)) * (x - (mu1 - mu2)), sigma_sq * sigma_sq); + hessian(p1, p2) -= w / sigma_sq; + hessian(p2, p1) -= w / sigma_sq; -#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(p1, p1) += w / sigma_sq; + hessian(p2, p2) += w / sigma_sq; + } +} - for (int j = 0; j < num_players; ++j) { - printf("%f ", fim[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"); } } @@ -197,9 +276,10 @@ int main(int argc, char **argv) for ( ;; ) { char pl1[256], pl2[256]; int score1, score2; + float weight; - if (scanf("%s %s %d %d", pl1, pl2, &score1, &score2) != 4) { - fprintf(stderr, "Read %d matches.\n", num_matches); + if (scanf("%s %s %d %d %f", pl1, pl2, &score1, &score2, &weight) != 5) { + //fprintf(stderr, "Read %d matches.\n", num_matches); break; } @@ -215,52 +295,67 @@ 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 < 100; ++j) { + for (int j = 0; j < 1000; ++j) { float old_mu[MAX_PLAYERS]; - float old_sigma[MAX_PLAYERS]; - memcpy(old_mu, mu, sizeof(float) * MAX_PLAYERS); - memcpy(old_sigma, sigma, sizeof(float) * MAX_PLAYERS); + float old_global_sigma = global_sigma; + float old_prior_sigma = prior_sigma; + memcpy(old_mu, mu, sizeof(mu)); 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, 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); } */ - bool any_difference = false; + + float sumdiff = 0.0f; for (int i = 0; i < num_players; ++i) { - if (fabs(mu[i] - old_mu[i]) > EPSILON || - fabs(sigma[i] - old_sigma[i]) > EPSILON) { - any_difference = true; - break; - } + sumdiff += (mu[i] - old_mu[i]) * (mu[i] - old_mu[i]); } - if (!any_difference) { - fprintf(stderr, "Converged after %d iterations. Stopping.\n", j); + 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("aux_param num_iterations %d\n", j + 1); break; } } - dump_scores(players, mu, sigma, num_players); -// construct_fim(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 }