X-Git-Url: https://git.sesse.net/?a=blobdiff_plain;f=bayeswf.cpp;h=487774ea394cb993c7b5c4499cd64f3fb6c30d55;hb=1e27ce64b3bd8574c621f1b7d47f94593d99318b;hp=1a4893050f4d3520479785f6ccec9cfa0819bad0;hpb=e6ad3fff4d7f562c9c68948b340777c90d4867c0;p=wloh diff --git a/bayeswf.cpp b/bayeswf.cpp index 1a48930..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,31 +33,37 @@ 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) { -#if 1 +#if 0 for (int i = 0; i < num_players; ++i) { - fprintf(stderr, "%s=[%5.1f, %4.1f] ", players[i].c_str(), mu[i], sigma[i]); + printf("%s=[%5.1f, %4.1f] ", players[i].c_str(), mu[i], sigma[i]); } - fprintf(stderr, "\n"); + printf("\n"); +#elif 0 + for (int i = 0; i < num_players; ++i) { + printf("%5.1f ", mu[i]); + } + printf("\n"); #else for (int i = 0; i < num_players; ++i) { - fprintf(stderr, "%5.1f ", mu[i]); + printf("%f %s\n", mu[i], players[i].c_str()); } - 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 ] + * 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) { @@ -62,10 +72,19 @@ 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]; - float inv_sigma_c2 = 1.0f / (sigma1 * sigma1 + sigma2 * sigma2); + float inv_sigma_c2 = matches[i].weight / (sigma1 * sigma1 + sigma2 * sigma2); float x = matches[i].margin; // / 70.0f; nom += (mu[matches[i].other_player] + x) * inv_sigma_c2; @@ -74,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 @@ -106,56 +143,126 @@ void update_sigma(float *mu, float *sigma, int player_num, const vector & sigma[player_num] = sqrt(sum / matches.size()); } -void renormalize(float *mu, float *sigma, int num_players) +/* + * 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_global_sigma(float *mu, float *sigma, int num_players) { - float avg = 0.0f; + float nom = 0.0f, denom = 0.0f; + 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. for (int i = 0; i < num_players; ++i) { - avg += mu[i]; + sigma[i] = best_sigma; } - float corr = 1500.0f - avg / 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 nom = 0.0f, denom = 0.0f; for (int i = 0; i < num_players; ++i) { - mu[i] += corr; + 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 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) { + 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; +} + /* - * 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: * - * The Hessian matrix is generally zero and thus not as interesting. + * M_ij = D_i D_j (- logL) = -w / sigma² for i != j + * w / sigma² for i == j + * + * 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) +double hessian[MAX_PLAYERS][MAX_PLAYERS]; +void construct_hessian(const float *mu, const float *sigma, int num_players) { - float fim[MAX_PLAYERS][MAX_PLAYERS]; - memset(fim, 0, sizeof(fim)); + memset(hessian, 0, sizeof(hessian)); - for (int i = 0; i < num_players; ++i) { - float mu1 = mu[i]; - float 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; - 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 sigma1 = sigma[m.player]; + double sigma2 = sigma[m.other_player]; - fprintf(stderr, "exp_diff_sq=%f sigma_sq=%f\n", (x - (mu1 - mu2)) * (x - (mu1 - mu2)), sigma_sq * sigma_sq); + double sigma_sq = sigma1 * sigma1 + sigma2 * sigma2; + float w = m.weight; -#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][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) { for (int j = 0; j < num_players; ++j) { - printf("%f ", fim[i][j]); + printf("%.12f ", hessian[i][j]); } printf("\n"); } @@ -192,9 +299,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; } @@ -210,54 +318,70 @@ 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); - - dump_scores(players, 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_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); - dump_scores(players, 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); } */ - 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]); + sumdiff += (sigma[i] - old_sigma[i]) * (sigma[i] - old_sigma[i]); } - if (!any_difference) { - fprintf(stderr, "Converged after %d iterations. Stopping.\n", j); + sumdiff += (prior_sigma - old_prior_sigma) * (prior_sigma - old_prior_sigma); + if (sumdiff < EPSILON) { + //fprintf(stderr, "Converged after %d iterations. Stopping.\n", j); + printf("%d -1\n", j + 1); break; } } -// construct_fim(mu, sigma, num_players); +#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)); + 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 }