--- /dev/null
+#include <stdio.h>
+#include <math.h>
+#include <string.h>
+#include <stdlib.h>
+
+#include <map>
+#include <vector>
+#include <string>
+#include <algorithm>
+
+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<int, vector<match> > matches_for_player;
+
+void dump_scores(const vector<string> &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<match> &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<match> &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<string> players;
+ map<string, int> 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);
+}
projects / wloh / commitdiff