13 #define MAX_PLAYERS 4096
15 float mu[MAX_PLAYERS];
16 float sigma[MAX_PLAYERS];
21 * L(mu_vec, sigma_vec, matches) = product[ L(mu_A, sigma_A, mu_B, sigma_B, score_AB - score_BA) ]
22 * log-likelihood = sum[ log( L(mu_A, sigma_A, mu_B, sigma_B, score_AB - score_BA) ) ]
24 * L(mu1, sigma1, mu2, sigma2, score2 - score1) = sigmoid(mu2 - mu1, sqrt(sigma1² + sigma2²), (score2 - score1))
26 * pdf := 1/(sigma * sqrt(2*Pi)) * exp(-(x - mu)^2 / (2 * sigma^2));
27 * pdfs := subs({ mu = mu1 - mu2, sigma = sqrt(sigma1^2 + sigma2^2) }, pdf);
28 * diff(log(pdfs), mu1);
36 map<int, vector<match> > matches_for_player;
38 void dump_scores(const vector<string> &players, const float *mu, const float *sigma, int num_players)
41 for (int i = 0; i < num_players; ++i) {
42 printf("%s=[%5.1f, %4.1f] ", players[i].c_str(), mu[i], sigma[i]);
46 for (int i = 0; i < num_players; ++i) {
47 printf("%5.1f ", mu[i]);
51 for (int i = 0; i < num_players; ++i) {
52 printf("%5.1f %s\n", mu[i], players[i].c_str());
59 * diff(logL, mu1) = -w * (mu1 - mu2 - x) / sigma_c^2
60 * maximizer for mu1 is given by: sum_i[ (w_i/sigma_c_i)^2 (mu1 - mu2_i - x_i) ] = 0
61 * sum_i[ (w_i/sigma_c_i)^2 mu1 ] = sum_i [ (w_i/sigma_c_i)^2 ( mu2_i + x_i ) ]
62 * mu1 = sum_i [ (w_i/sigma_c_i)^2 ( mu2_i + x_i ) ] / sum_i[ (w_i/sigma_c_i)^2 ]
64 void update_mu(float *mu, float *sigma, int player_num, const vector<match> &matches)
66 if (matches.empty()) {
70 float nom = 0.0f, denom = 0.0f;
71 for (unsigned i = 0; i < matches.size(); ++i) {
72 float sigma1 = sigma[player_num];
73 float sigma2 = sigma[matches[i].other_player];
74 float inv_sigma_c2 = matches[i].weight / (sigma1 * sigma1 + sigma2 * sigma2);
75 float x = matches[i].margin; // / 70.0f;
77 nom += (mu[matches[i].other_player] + x) * inv_sigma_c2;
78 denom += inv_sigma_c2;
80 mu[player_num] = nom / denom;
84 * diff(logL, sigma1) = sigma1 (-sigma1² - sigma2² + (x - mu)²) / sigma_c²
85 * maximizer for sigma1 is given by: sum_i[ (1/sigma_c_i)² sigma1 ((x - mu)² - (sigma1² + sigma2²) ] = 0
86 * sum_i[ (x - mu)² - sigma1² - sigma2² ] = 0 |: sigma1 != 0, sigma2 != 0
87 * sum_i[ (x - mu)² - sigma2² ] = sum[ sigma1² ]
88 * sigma1 = sqrt( sum_i[ (x - mu)² - sigma2² ] / N )
90 void update_sigma(float *mu, float *sigma, int player_num, const vector<match> &matches)
92 if (matches.size() < 2) {
97 for (unsigned i = 0; i < matches.size(); ++i) {
98 float mu1 = mu[player_num];
99 float mu2 = mu[matches[i].other_player];
100 float mu = mu1 - mu2;
101 float sigma2 = sigma[matches[i].other_player];
102 float x = matches[i].margin;
104 //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);
105 sum += (x - mu) * (x - mu) - sigma2 * sigma2;
111 //fprintf(stderr, "sum=%f\n", sum);
112 sigma[player_num] = sqrt(sum / matches.size());
116 * diff(logL, sigma) = w ( (x - mu)² - sigma² ) / sigma³
117 * maximizer for sigma is given by: sum_i[ (w_i/sigma)³ ((x - mu)² - sigma²) ] = 0
118 * sum_i[ w_i ( (x - mu)² - sigma² ) ] = 0 |: sigma != 0
119 * sum_i[ w_i (x - mu)² ] = sum[ w_i sigma² ]
120 * sigma = sqrt( sum_i[ w_i (x - mu)² ] / sum[w_i] )
122 void update_global_sigma(float *mu, float *sigma, int num_players)
124 float nom = 0.0f, denom = 0.0f;
125 for (int i = 0; i < num_players; ++i) {
126 for (unsigned j = 0; j < matches_for_player[i].size(); ++j) {
127 const match& m = matches_for_player[i][j];
129 // Only count each match once.
130 if (m.other_player <= i) {
135 float mu2 = mu[m.other_player];
136 float mu = mu1 - mu2;
140 nom += w * ((x - mu) * (x - mu));
145 float best_sigma = sqrt(nom / denom) / sqrt(2.0f); // Divide evenly between the two players.
146 for (int i = 0; i < num_players; ++i) {
147 sigma[i] = best_sigma;
151 void renormalize(float *mu, float *sigma, int num_players)
154 for (int i = 0; i < num_players; ++i) {
157 float corr = 1500.0f - avg / num_players;
158 for (int i = 0; i < num_players; ++i) {
164 * Compute Fisher information matrix of the log-likelihood, evaluated at the MLE,
166 * ie. M_ij = E[ (D_i logL) (D_j logL) ] = - sum( ( x - (mu_1 - mu_2) )² / sigma_c⁴ ) for i != j
167 * = - sum( 1 / sigma_c² ) for i == j
169 * The Hessian matrix is generally zero and thus not as interesting.
171 void construct_fim(const float *mu, const float *sigma, int num_players)
173 float fim[MAX_PLAYERS][MAX_PLAYERS];
174 memset(fim, 0, sizeof(fim));
176 for (int i = 0; i < num_players; ++i) {
178 float sigma1 = sigma[i];
180 for (unsigned k = 0; k < matches_for_player[i].size(); ++k) {
181 int j = matches_for_player[i][k].other_player;
183 float sigma2 = sigma[j];
185 float x = matches_for_player[i][k].margin;
186 float sigma_sq = sqrt(sigma1 * sigma1 + sigma2 * sigma2);
188 fprintf(stderr, "exp_diff_sq=%f sigma_sq=%f\n", (x - (mu1 - mu2)) * (x - (mu1 - mu2)), sigma_sq * sigma_sq);
191 fim[i][i] += (x - (mu1 - mu2)) * (x - (mu1 - mu2)) / (sigma_sq * sigma_sq);
192 fim[i][j] -= (x - (mu1 - mu2)) * (x - (mu1 - mu2)) / (sigma_sq * sigma_sq);
194 fim[i][i] -= 1.0f / sigma_sq;
195 fim[i][j] += 1.0f / sigma_sq;
199 for (int j = 0; j < num_players; ++j) {
200 printf("%f ", fim[i][j]);
206 int main(int argc, char **argv)
209 if (scanf("%d", &num_players) != 1) {
210 fprintf(stderr, "Could't read number of players\n");
214 if (num_players > MAX_PLAYERS) {
215 fprintf(stderr, "Max %d players supported\n", MAX_PLAYERS);
219 vector<string> players;
220 map<string, int> player_map;
222 for (int i = 0; i < num_players; ++i) {
224 if (scanf("%s", buf) != 1) {
225 fprintf(stderr, "Couldn't read player %d\n", i);
229 players.push_back(buf);
235 char pl1[256], pl2[256];
239 if (scanf("%s %s %d %d %f", pl1, pl2, &score1, &score2, &weight) != 5) {
240 fprintf(stderr, "Read %d matches.\n", num_matches);
246 if (player_map.count(pl1) == 0) {
247 fprintf(stderr, "Unknown player '%s'\n", pl1);
250 if (player_map.count(pl2) == 0) {
251 fprintf(stderr, "Unknown player '%s'\n", pl2);
256 m1.other_player = player_map[pl2];
257 m1.margin = score1 - score2;
259 matches_for_player[player_map[pl1]].push_back(m1);
262 m2.other_player = player_map[pl1];
263 m2.margin = score2 - score1;
265 matches_for_player[player_map[pl2]].push_back(m2);
268 float mu[MAX_PLAYERS];
269 float sigma[MAX_PLAYERS];
271 for (int i = 0; i < num_players; ++i) {
273 sigma[i] = 70.0f / sqrt(2.0f);
275 renormalize(mu, sigma, num_players);
277 for (int j = 0; j < 1000; ++j) {
278 float old_mu[MAX_PLAYERS];
279 float old_sigma[MAX_PLAYERS];
280 memcpy(old_mu, mu, sizeof(mu));
281 memcpy(old_sigma, sigma, sizeof(sigma));
282 for (int i = 0; i < num_players; ++i) {
283 update_mu(mu, sigma, i, matches_for_player[i]);
284 renormalize(mu, sigma, num_players);
286 update_global_sigma(mu, sigma, num_players);
287 /* for (int i = 0; i < num_players; ++i) {
288 update_sigma(mu, sigma, i, matches_for_player[i]);
289 dump_scores(players, mu, sigma, num_players);
292 float sumdiff = 0.0f;
293 for (int i = 0; i < num_players; ++i) {
294 sumdiff += (mu[i] - old_mu[i]) * (mu[i] - old_mu[i]);
295 sumdiff += (sigma[i] - old_sigma[i]) * (sigma[i] - old_sigma[i]);
297 if (sumdiff < EPSILON) {
298 fprintf(stderr, "Converged after %d iterations. Stopping.\n", j);
302 dump_scores(players, mu, sigma, num_players);
303 fprintf(stderr, "Optimal sigma: %f (two-player: %f)\n", sigma[0], sigma[0] * sqrt(2.0f));
305 // construct_fim(mu, sigma, num_players);