using namespace std;
+#define PRIOR_MU 1500
#define MAX_PLAYERS 4096
float mu[MAX_PLAYERS];
float sigma[MAX_PLAYERS];
+float prior_sigma = 70.0f;
#define EPSILON 1e-3
struct match {
int other_player;
int margin;
+ float weight;
};
map<int, vector<match> > matches_for_player;
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
}
/*
- * 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<match> &matches)
{
}
float nom = 0.0f, denom = 0.0f;
+
+ // Prior.
+ {
+ float inv_sigma2 = 1.0f / (prior_sigma * prior_sigma);
+ nom += PRIOR_MU * inv_sigma2;
+ denom += 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;
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 (int i = 0; i < num_players; ++i) {
- avg += mu[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;
+ }
+ }
+
+ 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;
}
- 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;
+ 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 w = m.weight;
+ nom += w * ((mu1 - PRIOR_MU) * (mu1 - PRIOR_MU));
+ denom += w * 1.0f;
+ }
}
+
+ prior_sigma = sqrt(nom / denom);
}
/*
- * 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];
+ double 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);
+ double sigma2 = sigma[j];
+ double sigma_sq = sigma1 * sigma1 + sigma2 * sigma2;
- fprintf(stderr, "exp_diff_sq=%f sigma_sq=%f\n", (x - (mu1 - mu2)) * (x - (mu1 - mu2)), sigma_sq * sigma_sq);
+ float w = matches_for_player[i][k].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[i][j] -= w / sigma_sq;
+ hessian[i][i] += 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");
}
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;
}
match m1;
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.other_player = player_map[pl1];
m2.margin = score2 - score1;
+ m2.weight = weight;
matches_for_player[player_map[pl2]].push_back(m2);
}
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) {
+ 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;
}
}
+ 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);
-// construct_fim(mu, sigma, num_players);
+// construct_hessian(mu, sigma, num_players);
}