Support the logistic distribution in addition to the normal distribution. Not enabled...

diff --git a/foosrank.cpp b/foosrank.cpp
index dd1f905f9cd92110d45b431eea8b2ff6545915e9..99f88b0d5909699630d382785bb59385f5757371 100644 (file)
--- a/foosrank.cpp
+++ b/foosrank.cpp
@@ -8,12 +8,19 @@
 #include <complex>
 #include <fftw3.h>
 
+#define USE_LOGISTIC_DISTRIBUTION 0
+
 // step sizes
 static const double int_step_size = 75.0;
 
 // rating constant (see below)
 static const double rating_constant = 455.0;
 
+#if USE_LOGISTIC_DISTRIBUTION
+// constant used in the logistic pdf
+static const double l_const = M_PI / (2.0 * sqrt(3.0));
+#endif
+
 using namespace std;
 
 static double prob_score(int k, int a, double rd);
@@ -21,6 +28,14 @@ static double prob_score_real(int k, int a, double binomial, double rd_norm);
 static double prodai(int k, int a);
 static double fac(int x);
 
+#if USE_LOGISTIC_DISTRIBUTION
+// sech²(x)
+static double sech2(double x)
+{
+       double c = cosh(x);
+       return 1.0 / (c*c);
+}
+#endif
 
 // probability of match ending k-a (k>a) when winnerR - loserR = RD
 //
@@ -122,11 +137,23 @@ static void compute_opponent_rating_pdf(int k, int a, double mu2, double sigma2,
                func1[i].real() = func1[i].imag() = func2[i].real() = func2[i].imag() = 0.0;
        }
 
+#if USE_LOGISTIC_DISTRIBUTION
+       double invsigma2 = 1.0 / sigma2;
+#else
        double invsq2sigma2 = 1.0 / (sqrt(2.0) * sigma2);
+#endif
        for (int i = 0; i < sz; ++i) {
                double x1 = 0.0 + h*i;
+
+               // opponent's pdf
+#if USE_LOGISTIC_DISTRIBUTION
+               double z = (x1 - mu2) * invsigma2;
+               double ch = cosh(l_const * z);
+               func1[i].real() = 1.0 / (ch * ch);
+#else
                double z = (x1 - mu2) * invsq2sigma2;
                func1[i].real() = exp(-z*z);
+#endif
 
                double x2 = -3000.0 + h*i;
                func2[(i - sz/2 + sz*2)%(sz*2)].real() = prob_score_real(k, a, binomial_precompute, x2*winfac);
@@ -286,17 +313,30 @@ static void least_squares(vector<pair<double, double> > &curve, double mu1, doub
                for (unsigned i = 0; i < curve.size(); ++i) {
                        double x = curve[i].first;
 
+#if USE_LOGISTIC_DISTRIBUTION
+                       // df/dA(x_i)
+                       matA[i + 0 * curve.size()] = sech2(l_const * (x-mu)/sigma);
+
+                       // df/dµ(x_i)
+                       matA[i + 1 * curve.size()] = 2.0 * l_const * A * matA[i + 0 * curve.size()]
+                               * tanh(l_const * (x-mu)/sigma) / sigma;
+
+                       // df/dσ(x_i)
+                       matA[i + 2 * curve.size()] = 
+                               matA[i + 1 * curve.size()] * (x-mu)/sigma;
+#else
                        // df/dA(x_i)
                        matA[i + 0 * curve.size()] = 
                                exp(-(x-mu)*(x-mu)/(2.0*sigma*sigma));
 
                        // df/dµ(x_i)
-                       matA[i + 1 * curve.size()] = 
+                       matA[i + 1 * curve.size()] =
                                A * (x-mu)/(sigma*sigma) * matA[i + 0 * curve.size()];
 
                        // df/dσ(x_i)
                        matA[i + 2 * curve.size()] = 
                                matA[i + 1 * curve.size()] * (x-mu)/sigma;
+#endif
                }
 
                // find dβ
@@ -304,7 +344,11 @@ static void least_squares(vector<pair<double, double> > &curve, double mu1, doub
                        double x = curve[i].first;
                        double y = curve[i].second;
 
+#if USE_LOGISTIC_DISTRIBUTION
+                       dbeta[i] = y - A * sech2(l_const * (x-mu)/sigma);
+#else
                        dbeta[i] = y - A * exp(- (x-mu)*(x-mu)/(2.0*sigma*sigma));
+#endif
                }
 
                // compute a and b
@@ -340,9 +384,16 @@ static void compute_new_rating(double mu1, double sigma1, double mu2, double sig
        // multiply in the gaussian
        for (unsigned i = 0; i < curve.size(); ++i) {
                double r1 = curve[i].first;
+
+               // my pdf
                double z = (r1 - mu1) / sigma1;
+#if USE_LOGISTIC_DISTRIBUTION
+               double ch = cosh(l_const * z);
+               curve[i].second /= (ch * ch);
+#else
                double gaussian = exp(-(z*z/2.0));
                curve[i].second *= gaussian;
+#endif
        }
 
        double mu_est, sigma_est;
@@ -375,18 +426,32 @@ static void compute_new_double_rating(double mu1, double sigma1, double mu2, dou
                double r1 = i * h;
 
                // iterate over r2
+#if USE_LOGISTIC_DISTRIBUTION
+               double invsigma2 = 1.0 / sigma2;
+#else
                double invsq2sigma2 = 1.0 / (sqrt(2.0) * sigma2);
+#endif
                for (unsigned j = 0; j < curve.size(); ++j) {
                        double r1plusr2 = curve[j].first;
                        double r2 = r1plusr2 - r1;
 
+#if USE_LOGISTIC_DISTRIBUTION
+                       double z = (r2 - mu2) * invsigma2;
+                       double gaussian = sech2(l_const * z);
+#else  
                        double z = (r2 - mu2) * invsq2sigma2;
                        double gaussian = exp(-z*z);
+#endif
                        sum += curve[j].second * gaussian;
                }
 
+#if USE_LOGISTIC_DISTRIBUTION
+               double z = (r1 - mu1) / sigma1;
+               double gaussian = sech2(l_const * z);
+#else
                double z = (r1 - mu1) / sigma1;
                double gaussian = exp(-(z*z/2.0));
+#endif
                newcurve.push_back(make_pair(r1, gaussian * sum));
        }