X-Git-Url: https://git.sesse.net/?a=blobdiff_plain;f=bayeswf.cpp;h=487774ea394cb993c7b5c4499cd64f3fb6c30d55;hb=1e27ce64b3bd8574c621f1b7d47f94593d99318b;hp=344300155ad7d68e4d964a4682b006758beba8f4;hpb=7e0f1e557809609e90b21018fd45994322b9bc26;p=wloh diff --git a/bayeswf.cpp b/bayeswf.cpp index 3443001..487774e 100644 --- a/bayeswf.cpp +++ b/bayeswf.cpp @@ -10,10 +10,14 @@ using namespace std; +#define PRIOR_MU 1500 +#define PRIOR_WEIGHT 1.0 #define MAX_PLAYERS 4096 +#define DUMP_RAW 0 float mu[MAX_PLAYERS]; float sigma[MAX_PLAYERS]; +float prior_sigma = 70.0f; #define EPSILON 1e-3 @@ -29,11 +33,12 @@ 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) { @@ -49,9 +54,8 @@ 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 %s\n", mu[i], players[i].c_str()); } - printf("\n"); #endif } @@ -68,6 +72,15 @@ void update_mu(float *mu, float *sigma, int player_num, const vector &mat } 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]; @@ -80,6 +93,24 @@ 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) +{ + 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 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 @@ -122,24 +153,17 @@ void update_sigma(float *mu, float *sigma, int player_num, const vector & void update_global_sigma(float *mu, float *sigma, 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 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; } float best_sigma = sqrt(nom / denom) / sqrt(2.0f); // Divide evenly between the two players. @@ -148,16 +172,60 @@ void update_global_sigma(float *mu, float *sigma, int num_players) } } -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, float *sigma, 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; } - float corr = 1500.0f - avg / num_players; + + prior_sigma = sqrt(nom / denom); + if (!(prior_sigma > 40.0f)) { + prior_sigma = 40.0f; + } +} + +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 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 sigma1 = sigma[m.player]; + 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; + + total_logl += w * compute_logl((x - mu) / sigma); + } + + return total_logl; } /* @@ -173,20 +241,23 @@ void construct_hessian(const float *mu, const float *sigma, int num_players) { memset(hessian, 0, sizeof(hessian)); - for (int i = 0; i < num_players; ++i) { - double sigma1 = sigma[i]; + 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 sigma1 = sigma[m.player]; + double sigma2 = sigma[m.other_player]; - float w = matches_for_player[i][k].weight; + double sigma_sq = sigma1 * sigma1 + sigma2 * sigma2; + float w = m.weight; - hessian[i][j] -= w / sigma_sq; - hessian[i][i] += w / sigma_sq; - } + hessian[p1][p2] -= w / sigma_sq; + hessian[p2][p1] -= w / sigma_sq; + + hessian[p1][p1] += w / sigma_sq; + hessian[p2][p2] += w / sigma_sq; } for (int i = 0; i < num_players; ++i) { @@ -231,7 +302,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 +318,41 @@ 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; + mu[i] = PRIOR_MU; sigma[i] = 70.0f / sqrt(2.0f); } - renormalize(mu, sigma, num_players); for (int j = 0; j < 1000; ++j) { float old_mu[MAX_PLAYERS]; float old_sigma[MAX_PLAYERS]; + 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_global_sigma(mu, sigma, num_players); + update_prior_sigma(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); @@ -288,13 +363,25 @@ int main(int argc, char **argv) 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); if (sumdiff < EPSILON) { - fprintf(stderr, "Converged after %d iterations. Stopping.\n", j); + //fprintf(stderr, "Converged after %d iterations. Stopping.\n", j); + printf("%d -1\n", j + 1); break; } } + +#if DUMP_RAW + dump_raw(mu, sigma, num_players); +#else dump_scores(players, mu, sigma, num_players); - fprintf(stderr, "Optimal sigma: %f (two-player: %f)\n", sigma[0], sigma[0] * sqrt(2.0f)); + //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); // construct_hessian(mu, sigma, num_players); +#endif }