#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);
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
//
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);
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β
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
// 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;
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));
}