From: Steinar H. Gunderson Date: Fri, 16 Mar 2012 23:23:05 +0000 (+0100) Subject: Initial commit. X-Git-Url: https://git.sesse.net/?a=commitdiff_plain;h=466dba7f5659cb4e68c13b2928ddbcab9e54d4b3;p=wloh Initial commit. --- 466dba7f5659cb4e68c13b2928ddbcab9e54d4b3 diff --git a/bayeswf.cpp b/bayeswf.cpp new file mode 100644 index 0000000..f745629 --- /dev/null +++ b/bayeswf.cpp @@ -0,0 +1,259 @@ +#include +#include +#include +#include + +#include +#include +#include +#include + +using namespace std; + +#define MAX_PLAYERS 16 +#define EPSILON 1e-3 + +/* + * L(mu_vec, sigma_vec, matches) = product[ L(mu_A, sigma_A, mu_B, sigma_B, score_AB - score_BA) ] + * log-likelihood = sum[ log( L(mu_A, sigma_A, mu_B, sigma_B, score_AB - score_BA) ) ] + * + * L(mu1, sigma1, mu2, sigma2, score2 - score1) = sigmoid(mu2 - mu1, sqrt(sigma1² + sigma2²), (score2 - score1)) + * + * pdf := 1/(sigma * sqrt(2*Pi)) * exp(-(x - mu)^2 / (2 * sigma^2)); + * pdfs := subs({ mu = mu1 - mu2, sigma = sqrt(sigma1^2 + sigma2^2) }, pdf); + * diff(log(pdfs), mu1); + */ + +struct match { + int other_player; + int margin; +}; +map > matches_for_player; + +void dump_scores(const vector &players, const float *mu, const float *sigma, int num_players) +{ +#if 1 + for (int i = 0; i < num_players; ++i) { + fprintf(stderr, "%s=[%5.1f, %4.1f] ", players[i].c_str(), mu[i], sigma[i]); + } + fprintf(stderr, "\n"); +#else + for (int i = 0; i < num_players; ++i) { + fprintf(stderr, "%5.1f ", mu[i]); + } + fprintf(stderr, "\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 ] + */ +void update_mu(float *mu, float *sigma, int player_num, const vector &matches) +{ + if (matches.empty()) { + return; + } + + float nom = 0.0f, denom = 0.0f; + 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 x = matches[i].margin; // / 70.0f; + + nom += (mu[matches[i].other_player] + x) * inv_sigma_c2; + denom += inv_sigma_c2; + } + mu[player_num] = nom / denom; +} + +/* + * 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()); +} + +void renormalize(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; + for (int i = 0; i < num_players; ++i) { + mu[i] += corr; + } +} + +/* + * 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 + * + * The Hessian matrix is generally zero and thus not as interesting. + */ +void construct_fim(const float *mu, const float *sigma, int num_players) +{ + float fim[MAX_PLAYERS][MAX_PLAYERS]; + memset(fim, 0, sizeof(fim)); + + for (int i = 0; i < num_players; ++i) { + float mu1 = mu[i]; + float sigma1 = sigma[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]; + + float x = matches_for_player[i][k].margin; + float sigma_sq = sqrt(sigma1 * sigma1 + sigma2 * sigma2); + + fprintf(stderr, "exp_diff_sq=%f sigma_sq=%f\n", (x - (mu1 - mu2)) * (x - (mu1 - mu2)), sigma_sq * 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 + } + + for (int j = 0; j < num_players; ++j) { + printf("%f ", fim[i][j]); + } + printf("\n"); + } +} + +int main(int argc, char **argv) +{ + int num_players; + if (scanf("%d", &num_players) != 1) { + fprintf(stderr, "Could't read number of players\n"); + exit(1); + } + + if (num_players > MAX_PLAYERS) { + fprintf(stderr, "Max %d players supported\n", MAX_PLAYERS); + exit(1); + } + + vector players; + map player_map; + + for (int i = 0; i < num_players; ++i) { + char buf[256]; + if (scanf("%s", buf) != 1) { + fprintf(stderr, "Couldn't read player %d\n", i); + exit(1); + } + + players.push_back(buf); + player_map[buf] = i; + } + + int num_matches = 0; + for ( ;; ) { + char pl1[256], pl2[256]; + int score1, score2; + + if (scanf("%s %s %d %d", pl1, pl2, &score1, &score2) != 4) { + fprintf(stderr, "Read %d matches.\n", num_matches); + break; + } + + ++num_matches; + + if (player_map.count(pl1) == 0) { + fprintf(stderr, "Unknown player '%s'\n", pl1); + exit(1); + } + if (player_map.count(pl2) == 0) { + fprintf(stderr, "Unknown player '%s'\n", pl2); + exit(1); + } + + match m1; + m1.other_player = player_map[pl2]; + m1.margin = score1 - score2; + matches_for_player[player_map[pl1]].push_back(m1); + + match m2; + m2.other_player = player_map[pl1]; + m2.margin = score2 - score1; + matches_for_player[player_map[pl2]].push_back(m2); + } + + 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); + } + renormalize(mu, sigma, num_players); + + dump_scores(players, mu, sigma, num_players); + + for (int j = 0; j < 100; ++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); + for (int i = 0; i < num_players; ++i) { + update_mu(mu, sigma, i, matches_for_player[i]); + renormalize(mu, sigma, num_players); + dump_scores(players, 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); + } */ + bool any_difference = false; + 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; + } + } + if (!any_difference) { + fprintf(stderr, "Converged after %d iterations. Stopping.\n", j); + break; + } + } + +// construct_fim(mu, sigma, num_players); +}