6 #include <Eigen/Eigenvalues>
14 using namespace Eigen;
17 #define PRIOR_WEIGHT 1.0
18 #define MAX_PLAYERS 4096
21 float mu[MAX_PLAYERS];
22 float mu_stddev[MAX_PLAYERS];
23 float global_sigma = 70.0f;
24 float prior_sigma = 70.0f;
29 * L(mu_vec, sigma_vec, matches) = product[ L(mu_A, sigma_A, mu_B, sigma_B, score_AB - score_BA) ]
30 * log-likelihood = sum[ log( L(mu_A, sigma_A, mu_B, sigma_B, score_AB - score_BA) ) ]
32 * L(mu1, sigma1, mu2, sigma2, score2 - score1) = sigmoid(mu2 - mu1, sqrt(sigma1² + sigma2²), (score2 - score1))
34 * pdf := 1/(sigma * sqrt(2*Pi)) * exp(-(x - mu)^2 / (2 * sigma^2));
35 * pdfs := subs({ mu = mu1 - mu2, sigma = sqrt(sigma1^2 + sigma2^2) }, pdf);
36 * diff(log(pdfs), mu1);
40 int player, other_player;
44 map<int, vector<match> > matches_for_player;
45 vector<match> all_matches;
47 void dump_scores(const vector<string> &players, const float *mu, const float *mu_stddev, int num_players)
50 for (int i = 0; i < num_players; ++i) {
51 printf("%s=[%5.1f, %4.1f] ", players[i].c_str(), mu[i], sigma[i]);
55 for (int i = 0; i < num_players; ++i) {
56 printf("%5.1f ", mu[i]);
60 for (int i = 0; i < num_players; ++i) {
61 printf("%f %f %s\n", mu[i], mu_stddev[i], players[i].c_str());
67 * diff(logL, mu1) = -w * (mu1 - mu2 - x) / sigma_c^2
68 * maximizer for mu1 is given by: sum_i[ (w_i/sigma_c_i)^2 (mu1 - mu2_i - x_i) ] = 0
69 * sum_i[ (w_i/sigma_c_i)^2 mu1 ] = sum_i [ (w_i/sigma_c_i)^2 ( mu2_i + x_i ) ]
70 * mu1 = sum_i [ (w_i/sigma_c_i)^2 ( mu2_i + x_i ) ] / sum_i[ (w_i/sigma_c_i)^2 ]
72 void update_mu(float *mu, int player_num, const vector<match> &matches)
74 if (matches.empty()) {
78 float nom = 0.0f, denom = 0.0f;
82 float inv_sigma2 = 1.0f / (prior_sigma * prior_sigma);
83 nom += PRIOR_WEIGHT * PRIOR_MU * inv_sigma2;
84 denom += PRIOR_WEIGHT * inv_sigma2;
88 for (unsigned i = 0; i < matches.size(); ++i) {
89 float inv_sigma_c2 = matches[i].weight / (global_sigma * global_sigma);
90 float x = matches[i].margin; // / 70.0f;
92 nom += (mu[matches[i].other_player] + x) * inv_sigma_c2;
93 denom += inv_sigma_c2;
95 mu[player_num] = nom / denom;
98 void dump_raw(const float *mu, int num_players)
100 for (unsigned i = 0; i < all_matches.size(); ++i) {
101 const match& m = all_matches[i];
103 float mu1 = mu[m.player];
104 float mu2 = mu[m.other_player];
105 float sigma = global_sigma;
106 float mu = mu1 - mu2;
110 printf("%f %f\n", (x - mu) / sigma, w);
115 * diff(logL, sigma) = w ( (x - mu)² - sigma² ) / sigma³
116 * maximizer for sigma is given by: sum_i[ (w_i/sigma)³ ((x - mu)² - sigma²) ] = 0
117 * sum_i[ w_i ( (x - mu)² - sigma² ) ] = 0 |: sigma != 0
118 * sum_i[ w_i (x - mu)² ] = sum[ w_i sigma² ]
119 * sigma = sqrt( sum_i[ w_i (x - mu)² ] / sum[w_i] )
121 void update_global_sigma(float *mu, int num_players)
123 float nom = 0.0f, denom = 0.0f;
124 for (unsigned i = 0; i < all_matches.size(); ++i) {
125 const match& m = all_matches[i];
127 float mu1 = mu[m.player];
128 float mu2 = mu[m.other_player];
129 float mu = mu1 - mu2;
133 nom += w * ((x - mu) * (x - mu));
137 global_sigma = sqrt(nom / denom);
141 * diff(priorlogL, sigma) = w ( (x - mu)² - sigma² ) / sigma³
142 * maximizer for sigma is given by: sum_i[ (w_i/sigma)³ ((x - mu)² - sigma²) ] = 0
143 * sum_i[ w_i ( (x - mu)² - sigma² ) ] = 0 |: sigma != 0
144 * sum_i[ w_i (x - mu)² ] = sum[ w_i sigma² ]
145 * sigma = sqrt( sum_i[ w_i (x - mu)² ] / sum[w_i] )
147 void update_prior_sigma(float *mu, int num_players)
149 float nom = 0.0f, denom = 0.0f;
150 for (int i = 0; i < num_players; ++i) {
153 nom += ((mu1 - PRIOR_MU) * (mu1 - PRIOR_MU));
157 prior_sigma = sqrt(nom / denom);
158 if (!(prior_sigma > 40.0f)) {
163 float compute_logl(float z)
165 return -0.5 * (log(2.0f / M_PI) + z * z);
168 float compute_total_logl(float *mu, int num_players)
170 float total_logl = 0.0f;
173 for (int i = 0; i < num_players; ++i) {
174 total_logl += PRIOR_WEIGHT * compute_logl((mu[i] - PRIOR_MU) / prior_sigma);
178 for (unsigned i = 0; i < all_matches.size(); ++i) {
179 const match& m = all_matches[i];
181 float mu1 = mu[m.player];
182 float mu2 = mu[m.other_player];
183 float sigma = global_sigma;
184 float mu = mu1 - mu2;
188 total_logl += w * compute_logl((x - mu) / sigma);
195 * Compute Hessian matrix of the negative log-likelihood, ie. for each term in logL:
197 * M_ij = D_i D_j (- logL) = -w / sigma² for i != j
198 * w / sigma² for i == j
200 * Note that this does not depend on mu or the margin at all.
202 Matrix<float, Dynamic, Dynamic> hessian;
203 void construct_hessian(const float *mu, int num_players)
205 hessian = Matrix<float, Dynamic, Dynamic>(num_players, num_players);
208 for (int i = 0; i < num_players; ++i) {
209 hessian(i, i) += 1.0f / (prior_sigma * prior_sigma);
211 for (unsigned i = 0; i < all_matches.size(); ++i) {
212 const match &m = all_matches[i];
215 int p2 = m.other_player;
217 double sigma_sq = global_sigma * global_sigma;
220 hessian(p1, p2) -= w / sigma_sq;
221 hessian(p2, p1) -= w / sigma_sq;
223 hessian(p1, p1) += w / sigma_sq;
224 hessian(p2, p2) += w / sigma_sq;
228 // Compute uncertainty (stddev) of mu estimates, which is sqrt((H^-1)_ii),
229 // where H is the Hessian (see construct_hessian()).
230 void compute_mu_uncertainty(const float *mu, const vector<string> &players)
232 // FIXME: Use pseudoinverse if applicable.
233 Matrix<float, Dynamic, Dynamic> ih = hessian.inverse();
234 for (unsigned i = 0; i < players.size(); ++i) {
235 mu_stddev[i] = sqrt(ih(i, i));
238 for (unsigned i = 0; i < players.size(); ++i) {
239 for (unsigned j = 0; j < players.size(); ++j) {
240 printf("covariance %s %s %f\n",
248 int main(int argc, char **argv)
251 if (scanf("%d", &num_players) != 1) {
252 fprintf(stderr, "Could't read number of players\n");
256 if (num_players > MAX_PLAYERS) {
257 fprintf(stderr, "Max %d players supported\n", MAX_PLAYERS);
261 vector<string> players;
262 map<string, int> player_map;
264 for (int i = 0; i < num_players; ++i) {
266 if (scanf("%s", buf) != 1) {
267 fprintf(stderr, "Couldn't read player %d\n", i);
271 players.push_back(buf);
277 char pl1[256], pl2[256];
281 if (scanf("%s %s %d %d %f", pl1, pl2, &score1, &score2, &weight) != 5) {
282 //fprintf(stderr, "Read %d matches.\n", num_matches);
288 if (player_map.count(pl1) == 0) {
289 fprintf(stderr, "Unknown player '%s'\n", pl1);
292 if (player_map.count(pl2) == 0) {
293 fprintf(stderr, "Unknown player '%s'\n", pl2);
298 m1.player = player_map[pl1];
299 m1.other_player = player_map[pl2];
300 m1.margin = score1 - score2;
302 matches_for_player[player_map[pl1]].push_back(m1);
305 m2.player = player_map[pl2];
306 m2.other_player = player_map[pl1];
307 m2.margin = score2 - score1;
309 matches_for_player[player_map[pl2]].push_back(m2);
311 all_matches.push_back(m1);
314 float mu[MAX_PLAYERS];
316 for (int i = 0; i < num_players; ++i) {
320 for (int j = 0; j < 1000; ++j) {
321 float old_mu[MAX_PLAYERS];
322 float old_global_sigma = global_sigma;
323 float old_prior_sigma = prior_sigma;
324 memcpy(old_mu, mu, sizeof(mu));
325 for (int i = 0; i < num_players; ++i) {
326 update_mu(mu, i, matches_for_player[i]);
328 update_global_sigma(mu, num_players);
329 update_prior_sigma(mu, num_players);
330 /* for (int i = 0; i < num_players; ++i) {
331 update_sigma(mu, sigma, i, matches_for_player[i]);
332 dump_scores(players, mu, sigma, num_players);
335 float sumdiff = 0.0f;
336 for (int i = 0; i < num_players; ++i) {
337 sumdiff += (mu[i] - old_mu[i]) * (mu[i] - old_mu[i]);
339 sumdiff += (prior_sigma - old_prior_sigma) * (prior_sigma - old_prior_sigma);
340 sumdiff += (global_sigma - old_global_sigma) * (global_sigma - old_global_sigma);
341 if (sumdiff < EPSILON) {
342 //fprintf(stderr, "Converged after %d iterations. Stopping.\n", j);
343 printf("aux_param num_iterations %d\n", j + 1);
349 dump_raw(mu, num_players);
351 construct_hessian(mu, num_players);
352 compute_mu_uncertainty(mu, players);
353 dump_scores(players, mu, mu_stddev, num_players);
354 //fprintf(stderr, "Optimal sigma: %f (two-player: %f)\n", sigma[0], sigma[0] * sqrt(2.0f));
355 printf("aux_param score_stddev %f\n", global_sigma / sqrt(2.0f));
356 printf("aux_param rating_prior_stddev %f\n", prior_sigma);
358 float total_logl = compute_total_logl(mu, num_players);
359 printf("aux_param total_log_likelihood %f\n", total_logl);