Add division view.

[wloh] / bayeswf.cpp

diff --git a/bayeswf.cpp b/bayeswf.cpp
index 176d3459e02d56baefd151702045662665c8c925..f2047bbf69549d510522bb316b33cf9788db3b5b 100644 (file)
--- a/bayeswf.cpp
+++ b/bayeswf.cpp
@@ -11,11 +11,13 @@
 using namespace std;
 
 #define PRIOR_MU 1500
-#define PRIOR_SIGMA 100
+#define PRIOR_WEIGHT 1.0
 #define MAX_PLAYERS 4096
+#define DUMP_RAW 0
 
 float mu[MAX_PLAYERS];
 float sigma[MAX_PLAYERS];
+float prior_sigma = 70.0f;
 
 #define EPSILON 1e-3
 
@@ -72,9 +74,9 @@ void update_mu(float *mu, float *sigma, int player_num, const vector<match> &mat
 
        // Prior.
        {
-               float inv_sigma2 = 1.0f / (PRIOR_SIGMA * PRIOR_SIGMA);
-               nom += PRIOR_MU * inv_sigma2;
-               denom += inv_sigma2;
+               float inv_sigma2 = 1.0f / (prior_sigma * prior_sigma);
+               nom += PRIOR_WEIGHT * PRIOR_MU * inv_sigma2;
+               denom += PRIOR_WEIGHT * inv_sigma2;
        }
 
        // All matches.
@@ -90,6 +92,31 @@ void update_mu(float *mu, float *sigma, int player_num, const vector<match> &mat
        mu[player_num] = nom / denom;
 }
 
+void dump_raw(const float *mu, const float *sigma, int num_players)
+{
+       for (int i = 0; i < num_players; ++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 sigma1 = sigma[i];
+                       float sigma2 = sigma[m.other_player];
+                       float sigma = sqrt(sigma1 * sigma1 + sigma2 * sigma2);
+                       float mu = mu1 - mu2;
+                       float x = m.margin;
+                       float w = m.weight;
+
+                       printf("%f %f\n", (x - mu) / sigma, w);
+               }
+       }
+}
+
 /*
  * 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
@@ -158,16 +185,77 @@ void update_global_sigma(float *mu, float *sigma, int num_players)
        }
 }
 
-void renormalize(float *mu, float *sigma, int 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) {
+               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);
+       if (!(prior_sigma > 40.0f)) {
+               prior_sigma = 40.0f;
+       }
+}
+
+float compute_logl(float z)
 {
-       float avg = 0.0f;
+       return -0.5 * (log(2.0f / M_PI) + z * z);
+}
+
+float compute_total_logl(float *mu, float *sigma, int num_players)
+{
+       float total_logl = 0.0f;
+
+       // Prior.
        for (int i = 0; i < num_players; ++i) {
-               avg += mu[i];
+               total_logl += PRIOR_WEIGHT * compute_logl((mu[i] - PRIOR_MU) / prior_sigma);
        }
-       float corr = 1500.0f - avg / num_players;
+
+       // Matches.
        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 mu2 = mu[m.other_player];
+                       float sigma1 = sigma[i];
+                       float sigma2 = sigma[m.other_player];
+                       float sigma = sqrt(sigma1 * sigma1 + sigma2 * sigma2);
+                       float mu = mu1 - mu2;
+                       float x = m.margin;
+                       float w = m.weight;
+
+                       total_logl += w * compute_logl((x - mu) / sigma);
+               }
        }
+
+       return total_logl;
 }
 
 /*
@@ -273,21 +361,21 @@ int main(int argc, char **argv)
        float sigma[MAX_PLAYERS];
 
        for (int i = 0; i < num_players; ++i) {
-               mu[i] = 1500.0f;
+               mu[i] = PRIOR_MU;
                sigma[i] = 70.0f / sqrt(2.0f);
        }
-       renormalize(mu, sigma, num_players);
 
        for (int j = 0; j < 1000; ++j) {
                float old_mu[MAX_PLAYERS];
                float old_sigma[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);
                }
                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);
@@ -298,15 +386,25 @@ int main(int argc, char **argv)
                        sumdiff += (mu[i] - old_mu[i]) * (mu[i] - old_mu[i]);
                        sumdiff += (sigma[i] - old_sigma[i]) * (sigma[i] - old_sigma[i]);
                }
+               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;
                }
        }
+
+#if DUMP_RAW
+       dump_raw(mu, sigma, num_players);
+#else
        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);
+
+       float total_logl = compute_total_logl(mu, sigma, num_players);
+       printf("%f -4\n", total_logl);
 
 //     construct_hessian(mu, sigma, num_players);
+#endif
 }