From 2e45e0e9f8d23adcf88fb89b7fa83a5a7d8d02ec Mon Sep 17 00:00:00 2001 From: "Steinar H. Gunderson" Date: Sun, 2 Dec 2007 12:43:15 +0100 Subject: [PATCH] Support the logistic distribution in addition to the normal distribution. Not enabled yet. --- foosrank.cpp | 67 +++++++++++++++++++++++++++++++++++++++++++++++++++- 1 file changed, 66 insertions(+), 1 deletion(-) diff --git a/foosrank.cpp b/foosrank.cpp index dd1f905..99f88b0 100644 --- a/foosrank.cpp +++ b/foosrank.cpp @@ -8,12 +8,19 @@ #include #include +#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 > &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 > &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)); } -- 2.39.2