14 #define PRIOR_WEIGHT 1.0
15 #define MAX_PLAYERS 4096
18 float mu[MAX_PLAYERS];
19 float sigma[MAX_PLAYERS];
20 float prior_sigma = 70.0f;
25 * L(mu_vec, sigma_vec, matches) = product[ L(mu_A, sigma_A, mu_B, sigma_B, score_AB - score_BA) ]
26 * log-likelihood = sum[ log( L(mu_A, sigma_A, mu_B, sigma_B, score_AB - score_BA) ) ]
28 * L(mu1, sigma1, mu2, sigma2, score2 - score1) = sigmoid(mu2 - mu1, sqrt(sigma1² + sigma2²), (score2 - score1))
30 * pdf := 1/(sigma * sqrt(2*Pi)) * exp(-(x - mu)^2 / (2 * sigma^2));
31 * pdfs := subs({ mu = mu1 - mu2, sigma = sqrt(sigma1^2 + sigma2^2) }, pdf);
32 * diff(log(pdfs), mu1);
40 map<int, vector<match> > matches_for_player;
42 void dump_scores(const vector<string> &players, const float *mu, const float *sigma, int num_players)
45 for (int i = 0; i < num_players; ++i) {
46 printf("%s=[%5.1f, %4.1f] ", players[i].c_str(), mu[i], sigma[i]);
50 for (int i = 0; i < num_players; ++i) {
51 printf("%5.1f ", mu[i]);
55 for (int i = 0; i < num_players; ++i) {
56 printf("%f %s\n", mu[i], players[i].c_str());
62 * diff(logL, mu1) = -w * (mu1 - mu2 - x) / sigma_c^2
63 * maximizer for mu1 is given by: sum_i[ (w_i/sigma_c_i)^2 (mu1 - mu2_i - x_i) ] = 0
64 * sum_i[ (w_i/sigma_c_i)^2 mu1 ] = sum_i [ (w_i/sigma_c_i)^2 ( mu2_i + x_i ) ]
65 * mu1 = sum_i [ (w_i/sigma_c_i)^2 ( mu2_i + x_i ) ] / sum_i[ (w_i/sigma_c_i)^2 ]
67 void update_mu(float *mu, float *sigma, int player_num, const vector<match> &matches)
69 if (matches.empty()) {
73 float nom = 0.0f, denom = 0.0f;
77 float inv_sigma2 = 1.0f / (prior_sigma * prior_sigma);
78 nom += PRIOR_WEIGHT * PRIOR_MU * inv_sigma2;
79 denom += PRIOR_WEIGHT * inv_sigma2;
83 for (unsigned i = 0; i < matches.size(); ++i) {
84 float sigma1 = sigma[player_num];
85 float sigma2 = sigma[matches[i].other_player];
86 float inv_sigma_c2 = matches[i].weight / (sigma1 * sigma1 + sigma2 * sigma2);
87 float x = matches[i].margin; // / 70.0f;
89 nom += (mu[matches[i].other_player] + x) * inv_sigma_c2;
90 denom += inv_sigma_c2;
92 mu[player_num] = nom / denom;
95 void dump_raw(const float *mu, const float *sigma, int num_players)
97 for (int i = 0; i < num_players; ++i) {
98 for (unsigned j = 0; j < matches_for_player[i].size(); ++j) {
99 const match& m = matches_for_player[i][j];
101 // Only count each match once.
102 if (m.other_player <= i) {
107 float mu2 = mu[m.other_player];
108 float sigma1 = sigma[i];
109 float sigma2 = sigma[m.other_player];
110 float sigma = sqrt(sigma1 * sigma1 + sigma2 * sigma2);
111 float mu = mu1 - mu2;
115 printf("%f %f\n", (x - mu) / sigma, w);
121 * diff(logL, sigma1) = sigma1 (-sigma1² - sigma2² + (x - mu)²) / sigma_c²
122 * maximizer for sigma1 is given by: sum_i[ (1/sigma_c_i)² sigma1 ((x - mu)² - (sigma1² + sigma2²) ] = 0
123 * sum_i[ (x - mu)² - sigma1² - sigma2² ] = 0 |: sigma1 != 0, sigma2 != 0
124 * sum_i[ (x - mu)² - sigma2² ] = sum[ sigma1² ]
125 * sigma1 = sqrt( sum_i[ (x - mu)² - sigma2² ] / N )
127 void update_sigma(float *mu, float *sigma, int player_num, const vector<match> &matches)
129 if (matches.size() < 2) {
134 for (unsigned i = 0; i < matches.size(); ++i) {
135 float mu1 = mu[player_num];
136 float mu2 = mu[matches[i].other_player];
137 float mu = mu1 - mu2;
138 float sigma2 = sigma[matches[i].other_player];
139 float x = matches[i].margin;
141 //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);
142 sum += (x - mu) * (x - mu) - sigma2 * sigma2;
148 //fprintf(stderr, "sum=%f\n", sum);
149 sigma[player_num] = sqrt(sum / matches.size());
153 * diff(logL, sigma) = w ( (x - mu)² - sigma² ) / sigma³
154 * maximizer for sigma is given by: sum_i[ (w_i/sigma)³ ((x - mu)² - sigma²) ] = 0
155 * sum_i[ w_i ( (x - mu)² - sigma² ) ] = 0 |: sigma != 0
156 * sum_i[ w_i (x - mu)² ] = sum[ w_i sigma² ]
157 * sigma = sqrt( sum_i[ w_i (x - mu)² ] / sum[w_i] )
159 void update_global_sigma(float *mu, float *sigma, int num_players)
161 float nom = 0.0f, denom = 0.0f;
162 for (int i = 0; i < num_players; ++i) {
163 for (unsigned j = 0; j < matches_for_player[i].size(); ++j) {
164 const match& m = matches_for_player[i][j];
166 // Only count each match once.
167 if (m.other_player <= i) {
172 float mu2 = mu[m.other_player];
173 float mu = mu1 - mu2;
177 nom += w * ((x - mu) * (x - mu));
182 float best_sigma = sqrt(nom / denom) / sqrt(2.0f); // Divide evenly between the two players.
183 for (int i = 0; i < num_players; ++i) {
184 sigma[i] = best_sigma;
189 * diff(priorlogL, sigma) = w ( (x - mu)² - sigma² ) / sigma³
190 * maximizer for sigma is given by: sum_i[ (w_i/sigma)³ ((x - mu)² - sigma²) ] = 0
191 * sum_i[ w_i ( (x - mu)² - sigma² ) ] = 0 |: sigma != 0
192 * sum_i[ w_i (x - mu)² ] = sum[ w_i sigma² ]
193 * sigma = sqrt( sum_i[ w_i (x - mu)² ] / sum[w_i] )
195 void update_prior_sigma(float *mu, float *sigma, int num_players)
197 float nom = 0.0f, denom = 0.0f;
198 for (int i = 0; i < num_players; ++i) {
199 for (unsigned j = 0; j < matches_for_player[i].size(); ++j) {
200 const match& m = matches_for_player[i][j];
202 // Only count each match once.
203 if (m.other_player <= i) {
210 nom += w * ((mu1 - PRIOR_MU) * (mu1 - PRIOR_MU));
215 prior_sigma = sqrt(nom / denom);
216 if (!(prior_sigma > 40.0f)) {
221 float compute_logl(float z)
223 return -0.5 * (log(2.0f / M_PI) + z * z);
226 float compute_total_logl(float *mu, float *sigma, int num_players)
228 float total_logl = 0.0f;
231 for (int i = 0; i < num_players; ++i) {
232 total_logl += PRIOR_WEIGHT * compute_logl((mu[i] - PRIOR_MU) / prior_sigma);
236 for (int i = 0; i < num_players; ++i) {
237 for (unsigned j = 0; j < matches_for_player[i].size(); ++j) {
238 const match& m = matches_for_player[i][j];
240 // Only count each match once.
241 if (m.other_player <= i) {
246 float mu2 = mu[m.other_player];
247 float sigma1 = sigma[i];
248 float sigma2 = sigma[m.other_player];
249 float sigma = sqrt(sigma1 * sigma1 + sigma2 * sigma2);
250 float mu = mu1 - mu2;
254 total_logl += w * compute_logl((x - mu) / sigma);
262 * Compute Hessian matrix of the negative log-likelihood, ie. for each term in logL:
264 * M_ij = D_i D_j (- logL) = -w / sigma² for i != j
265 * w / sigma² for i == j
267 * Note that this does not depend on mu or the margin at all.
269 double hessian[MAX_PLAYERS][MAX_PLAYERS];
270 void construct_hessian(const float *mu, const float *sigma, int num_players)
272 memset(hessian, 0, sizeof(hessian));
274 for (int i = 0; i < num_players; ++i) {
275 double sigma1 = sigma[i];
277 for (unsigned k = 0; k < matches_for_player[i].size(); ++k) {
278 int j = matches_for_player[i][k].other_player;
280 double sigma2 = sigma[j];
281 double sigma_sq = sigma1 * sigma1 + sigma2 * sigma2;
283 float w = matches_for_player[i][k].weight;
285 hessian[i][j] -= w / sigma_sq;
286 hessian[i][i] += w / sigma_sq;
290 for (int i = 0; i < num_players; ++i) {
291 for (int j = 0; j < num_players; ++j) {
292 printf("%.12f ", hessian[i][j]);
298 int main(int argc, char **argv)
301 if (scanf("%d", &num_players) != 1) {
302 fprintf(stderr, "Could't read number of players\n");
306 if (num_players > MAX_PLAYERS) {
307 fprintf(stderr, "Max %d players supported\n", MAX_PLAYERS);
311 vector<string> players;
312 map<string, int> player_map;
314 for (int i = 0; i < num_players; ++i) {
316 if (scanf("%s", buf) != 1) {
317 fprintf(stderr, "Couldn't read player %d\n", i);
321 players.push_back(buf);
327 char pl1[256], pl2[256];
331 if (scanf("%s %s %d %d %f", pl1, pl2, &score1, &score2, &weight) != 5) {
332 //fprintf(stderr, "Read %d matches.\n", num_matches);
338 if (player_map.count(pl1) == 0) {
339 fprintf(stderr, "Unknown player '%s'\n", pl1);
342 if (player_map.count(pl2) == 0) {
343 fprintf(stderr, "Unknown player '%s'\n", pl2);
348 m1.other_player = player_map[pl2];
349 m1.margin = score1 - score2;
351 matches_for_player[player_map[pl1]].push_back(m1);
354 m2.other_player = player_map[pl1];
355 m2.margin = score2 - score1;
357 matches_for_player[player_map[pl2]].push_back(m2);
360 float mu[MAX_PLAYERS];
361 float sigma[MAX_PLAYERS];
363 for (int i = 0; i < num_players; ++i) {
365 sigma[i] = 70.0f / sqrt(2.0f);
368 for (int j = 0; j < 1000; ++j) {
369 float old_mu[MAX_PLAYERS];
370 float old_sigma[MAX_PLAYERS];
371 float old_prior_sigma = prior_sigma;
372 memcpy(old_mu, mu, sizeof(mu));
373 memcpy(old_sigma, sigma, sizeof(sigma));
374 for (int i = 0; i < num_players; ++i) {
375 update_mu(mu, sigma, i, matches_for_player[i]);
377 update_global_sigma(mu, sigma, num_players);
378 update_prior_sigma(mu, sigma, num_players);
379 /* for (int i = 0; i < num_players; ++i) {
380 update_sigma(mu, sigma, i, matches_for_player[i]);
381 dump_scores(players, mu, sigma, num_players);
384 float sumdiff = 0.0f;
385 for (int i = 0; i < num_players; ++i) {
386 sumdiff += (mu[i] - old_mu[i]) * (mu[i] - old_mu[i]);
387 sumdiff += (sigma[i] - old_sigma[i]) * (sigma[i] - old_sigma[i]);
389 sumdiff += (prior_sigma - old_prior_sigma) * (prior_sigma - old_prior_sigma);
390 if (sumdiff < EPSILON) {
391 //fprintf(stderr, "Converged after %d iterations. Stopping.\n", j);
392 printf("%d -1\n", j + 1);
398 dump_raw(mu, sigma, num_players);
400 dump_scores(players, mu, sigma, num_players);
401 //fprintf(stderr, "Optimal sigma: %f (two-player: %f)\n", sigma[0], sigma[0] * sqrt(2.0f));
402 printf("%f -2\n", sigma[0]);
403 printf("%f -3\n", prior_sigma);
405 float total_logl = compute_total_logl(mu, sigma, num_players);
406 printf("%f -4\n", total_logl);
408 // construct_hessian(mu, sigma, num_players);