#include <vector>
#include <algorithm>
-// integration step size
-static const double step_size = 10.0;
+// step sizes
+static const double int_step_size = 50.0;
+static const double pdf_step_size = 10.0;
+
+// rating constant (see below)
+static const double rating_constant = 455.0;
using namespace std;
//
double prob_score(double a, double rd)
{
- return prob_score_real(a, prodai(a), rd/455.0);
+ return prob_score_real(a, prodai(a), rd/rating_constant);
}
// Same, but takes in Product(a+i, i=1..9) as an argument in addition to a. Faster
// Set the last parameter to 1.0 if player 1 won, or -1.0 if player 2 won.
// In the latter case, ProbScore will be given (r1-r2) instead of (r2-r1).
//
+static inline double evaluate_int_point(double a, double prodai_precompute, double r1, double mu2, double sigma2, double winfac, double x);
+
double opponent_rating_pdf(double a, double r1, double mu2, double sigma2, double winfac)
{
- double sum = 0.0;
double prodai_precompute = prodai(a);
- winfac /= 455.0;
- for (double r2 = 0.0; r2 < 3000.0; r2 += step_size) {
- double x = r2 + step_size*0.5;
- double probscore = prob_score_real(a, prodai_precompute, (r1 - x)*winfac);
- double z = (x - mu2)/sigma2;
- double gaussian = exp(-(z*z/2.0));
- sum += step_size * probscore * gaussian;
+ winfac /= rating_constant;
+
+ int n = int(3000.0 / int_step_size + 0.5);
+ double h = 3000.0 / double(n);
+ double sum = evaluate_int_point(a, prodai_precompute, r1, mu2, sigma2, winfac, 0.0);
+
+ for (int i = 1; i < n; i += 2) {
+ sum += 4.0 * evaluate_int_point(a, prodai_precompute, r1, mu2, sigma2, winfac, i * h);
}
- return sum;
+ for (int i = 2; i < n; i += 2) {
+ sum += 2.0 * evaluate_int_point(a, prodai_precompute, r1, mu2, sigma2, winfac, i * h);
+ }
+ sum += evaluate_int_point(a, prodai_precompute, r1, mu2, sigma2, winfac, 3000.0);
+
+ return (h/3.0) * sum;
+}
+
+static inline double evaluate_int_point(double a, double prodai_precompute, double r1, double mu2, double sigma2, double winfac, double x)
+{
+ double probscore = prob_score_real(a, prodai_precompute, (r1 - x)*winfac);
+ double z = (x - mu2)/sigma2;
+ double gaussian = exp(-(z*z/2.0));
+ return probscore * gaussian;
}
// normalize the curve so we know that A ~= 1
}
// Give an OK starting estimate for the least squares, by numerical integration
-// of x*f(x) and x^2 * f(x). Somehow seems to underestimate sigma, though.
+// of statistical moments.
void estimate_musigma(vector<pair<double, double> > &curve, double &mu_result, double &sigma_result)
{
- double mu = 0.0;
- double sigma = 0.0;
- double sum_area = 0.0;
-
- for (unsigned i = 1; i < curve.size(); ++i) {
- double x1 = curve[i].first;
- double x0 = curve[i-1].first;
- double y1 = curve[i].second;
- double y0 = curve[i-1].second;
- double xm = 0.5 * (x0 + x1);
- double ym = 0.5 * (y0 + y1);
- sum_area += (x1-x0) * ym;
- mu += (x1-x0) * xm * ym;
- sigma += (x1-x0) * xm * xm * ym;
+ double h = (curve.back().first - curve.front().first) / (curve.size() - 1);
+
+ double area = curve.front().second;
+ double ex = curve.front().first * curve.front().second;
+ double ex2 = curve.front().first * curve.front().first * curve.front().second;
+
+ for (unsigned i = 1; i < curve.size() - 1; i += 2) {
+ double x = curve[i].first;
+ double y = curve[i].second;
+ area += 4.0 * y;
+ ex += 4.0 * x * y;
+ ex2 += 4.0 * x * x * y;
+ }
+ for (unsigned i = 2; i < curve.size() - 1; i += 2) {
+ double x = curve[i].first;
+ double y = curve[i].second;
+ area += 2.0 * y;
+ ex += 2.0 * x * y;
+ ex2 += 2.0 * x * x * y;
}
+
+ area += curve.back().second;
+ ex += curve.back().first * curve.back().second;
+ ex2 += curve.back().first * curve.back().first * curve.back().second;
+
+ area = (h/3.0) * area;
+ ex = (h/3.0) * ex / area;
+ ex2 = (h/3.0) * ex2 / area;
- mu_result = mu / sum_area;
- sigma_result = sqrt(sigma) / sum_area;
+ mu_result = ex;
+ sigma_result = sqrt(ex2 - ex * ex);
}
// Find best fit of the data in curves to a Gaussian pdf, based on the
// given initial estimates. Works by nonlinear least squares, iterating
// until we're below a certain threshold.
+//
+// Note that the algorithm blows up quite hard if the initial estimate is
+// not good enough. Use estimate_musigma to get a reasonable starting
+// estimate.
void least_squares(vector<pair<double, double> > &curve, double mu1, double sigma1, double &mu_result, double &sigma_result)
{
double A = 1.0;
sigma_result = sigma;
}
-int main(int argc, char **argv)
+void compute_new_rating(double mu1, double sigma1, double mu2, double sigma2, int score1, int score2, double &mu, double &sigma)
{
- double mu1 = atof(argv[1]);
- double sigma1 = atof(argv[2]);
- double mu2 = atof(argv[3]);
- double sigma2 = atof(argv[4]);
- int score1 = atoi(argv[5]);
- int score2 = atoi(argv[6]);
vector<pair<double, double> > curve;
if (score1 == 10) {
- for (double r1 = 0.0; r1 < 3000.0; r1 += step_size) {
+ for (double r1 = 0.0; r1 < 3000.0; r1 += pdf_step_size) {
double z = (r1 - mu1) / sigma1;
double gaussian = exp(-(z*z/2.0));
curve.push_back(make_pair(r1, gaussian * opponent_rating_pdf(score2, r1, mu2, sigma2, 1.0)));
}
} else {
- for (double r1 = 0.0; r1 < 3000.0; r1 += step_size) {
+ for (double r1 = 0.0; r1 < 3000.0; r1 += pdf_step_size) {
double z = (r1 - mu1) / sigma1;
double gaussian = exp(-(z*z/2.0));
curve.push_back(make_pair(r1, gaussian * opponent_rating_pdf(score1, r1, mu2, sigma2, -1.0)));
}
}
- double mu_est, sigma_est, mu, sigma;
+ double mu_est, sigma_est;
normalize(curve);
estimate_musigma(curve, mu_est, sigma_est);
least_squares(curve, mu_est, sigma_est, mu, sigma);
+}
+
+int main(int argc, char **argv)
+{
+ double mu1 = atof(argv[1]);
+ double sigma1 = atof(argv[2]);
+ double mu2 = atof(argv[3]);
+ double sigma2 = atof(argv[4]);
+ int score1 = atoi(argv[5]);
+ int score2 = atoi(argv[6]);
+ double mu, sigma;
+ compute_new_rating(mu1, sigma1, mu2, sigma2, score1, score2, mu, sigma);
printf("%f %f\n", mu, sigma);
}
+
projects / foosball / blobdiff