// integration step size
static const double 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);
//
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
{
double sum = 0.0;
double prodai_precompute = prodai(a);
- winfac /= 455.0;
+ winfac /= rating_constant;
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);
}
// 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;
+ double ex = 0.0;
+ double ex2 = 0.0;
for (unsigned i = 1; i < curve.size(); ++i) {
double x1 = curve[i].first;
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;
+ ex += (x1-x0) * xm * ym;
+ ex2 += (x1-x0) * xm * xm * ym;
}
- mu_result = mu / sum_area;
- sigma_result = sqrt(sigma) / sum_area;
+ ex /= sum_area;
+ ex2 /= 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;