// 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_real(int k, int a, double binomial, double rd_norm);
// 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);
+ func1[i].real() = sech2(0.5 * z);
#else
double z = (x1 - mu2) * invsq2sigma2;
func1[i].real() = exp(-z*z);
#if USE_LOGISTIC_DISTRIBUTION
// df/dA(x_i)
- matA[i + 0 * curve.size()] = sech2(l_const * (x-mu)/sigma);
+ matA[i + 0 * curve.size()] = sech2(0.5 * (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;
+ matA[i + 1 * curve.size()] = A * matA[i + 0 * curve.size()]
+ * tanh(0.5 * (x-mu)/sigma) / sigma;
// df/dσ(x_i)
matA[i + 2 * curve.size()] =
double y = curve[i].second;
#if USE_LOGISTIC_DISTRIBUTION
- dbeta[i] = y - A * sech2(l_const * (x-mu)/sigma);
+ dbeta[i] = y - A * sech2(0.5 * (x-mu)/sigma);
#else
dbeta[i] = y - A * exp(- (x-mu)*(x-mu)/(2.0*sigma*sigma));
#endif
// my pdf
double z = (r1 - mu1) / sigma1;
#if USE_LOGISTIC_DISTRIBUTION
- double ch = cosh(l_const * z);
- curve[i].second /= (ch * ch);
+ curve[i].second *= sech2(0.5 * z);
#else
double gaussian = exp(-(z*z/2.0));
curve[i].second *= gaussian;
#endif
}
-
+
// Compute the overall probability of the given result, by integrating
// the entire resulting pdf. Note that since we're actually evaluating
// a double integral, we'll need to multiply by h² instead of h.
- // (For some reason, Simpson's rule gives markedly worse accuracy here.
- // Probably related to h² somehow?)
{
double h = (curve.back().first - curve.front().first) / (curve.size() - 1);
- double sum = 0.0;
- for (unsigned i = 0; i < curve.size(); ++i) {
- sum += curve[i].second;
+ double sum = curve.front().second;
+ for (unsigned i = 1; i < curve.size() - 1; i += 2) {
+ sum += 4.0 * curve[i].second;
}
-
- sum *= h * h;
+ for (unsigned i = 2; i < curve.size() - 1; i += 2) {
+ sum += 2.0 * curve[i].second;
+ }
+ sum += curve.back().second;
+ sum *= h * h / 3.0;
// FFT convolution multiplication factor (FFTW computes unnormalized
// transforms)
// pdf normalization factors
#if USE_LOGISTIC_DISTRIBUTION
- sum *= M_PI / (sigma1 * 4.0 * sqrt(3.0));
- sum *= M_PI / (sigma2 * 4.0 * sqrt(3.0));
+ sum /= (sigma1 * 4.0);
+ sum /= (sigma2 * 4.0);
#else
sum /= (sigma1 * sqrt(2.0 * M_PI));
sum /= (sigma2 * sqrt(2.0 * M_PI));
#if USE_LOGISTIC_DISTRIBUTION
double z = (r2 - mu2) * invsigma2;
- double gaussian = sech2(l_const * z);
+ double gaussian = sech2(0.5 * z);
#else
double z = (r2 - mu2) * invsq2sigma2;
double gaussian = exp(-z*z);
#if USE_LOGISTIC_DISTRIBUTION
double z = (r1 - mu1) / sigma1;
- double gaussian = sech2(l_const * z);
+ double gaussian = sech2(0.5 * z);
#else
double z = (r1 - mu1) / sigma1;
double gaussian = exp(-(z*z/2.0));
// a triple integral, we'll need to multiply by 4h³ (no idea where the
// 4 factor comes from, probably from the 0..6000 range somehow) instead
// of h.
- // (TODO: Evaluate Simpson's rule here, although it's probably even worse
- // than for the single case.)
{
double h = (newcurve.back().first - newcurve.front().first) / (newcurve.size() - 1);
- double sum = 0.0;
- for (unsigned i = 0; i < newcurve.size(); ++i) {
- sum += newcurve[i].second;
+ double sum = newcurve.front().second;
+ for (unsigned i = 1; i < newcurve.size() - 1; i += 2) {
+ sum += 4.0 * newcurve[i].second;
}
+ for (unsigned i = 2; i < newcurve.size() - 1; i += 2) {
+ sum += 2.0 * newcurve[i].second;
+ }
+ sum += newcurve.back().second;
- sum *= 4.0 * h * h * h;
+ sum *= 4.0 * h * h * h / 3.0;
// FFT convolution multiplication factor (FFTW computes unnormalized
// transforms)
// pdf normalization factors
#if USE_LOGISTIC_DISTRIBUTION
- sum *= M_PI / (sigma1 * 4.0 * sqrt(3.0));
- sum *= M_PI / (sigma2 * 4.0 * sqrt(3.0));
- sum *= M_PI / (sigma_t * 4.0 * sqrt(3.0));
+ sum /= (sigma1 * 4.0);
+ sum /= (sigma2 * 4.0);
+ sum /= (sigma_t * 4.0);
#else
sum /= (sigma1 * sqrt(2.0 * M_PI));
sum /= (sigma2 * sqrt(2.0 * M_PI));
int score2 = atoi(argv[10]);
double mu, sigma, probability;
compute_new_double_rating(mu1, sigma1, mu2, sigma2, mu3, sigma3, mu4, sigma4, score1, score2, mu, sigma, probability);
- if (score1 > score2) {
- printf("%f %f %f\n", mu, sigma, probability);
- } else {
- printf("%f %f %f\n", mu, sigma, probability);
- }
+ printf("%f %f %f\n", mu, sigma, probability);
} else if (argc > 8) {
double mu3 = atof(argv[5]);
double sigma3 = atof(argv[6]);
double mu4 = atof(argv[7]);
double sigma4 = atof(argv[8]);
int k = atoi(argv[9]);
-
+
// assess all possible scores
for (int i = 0; i < k; ++i) {
double newmu1_1, newmu1_2, newmu2_1, newmu2_2;
double mu, sigma, probability;
compute_new_rating(mu1, sigma1, mu2, sigma2, score1, score2, mu, sigma, probability);
- if (score1 > score2) {
- printf("%f %f %f\n", mu, sigma, probability);
- } else {
- printf("%f %f %f\n", mu, sigma, probability);
- }
+ printf("%f %f %f\n", mu, sigma, probability);
} else {
int k = atoi(argv[5]);
projects / foosball / blobdiff