X-Git-Url: https://git.sesse.net/?a=blobdiff_plain;f=foorank.cpp;h=f8ecf63872be9f5abd67af011f531ba595a901b6;hb=130255c3b613130f3799dfa17ece66b5c6d62e06;hp=096c185f3db0ad1ed8159e6f8ac70c3145c0c615;hpb=2d652fe32a6b7e994050be7faad1e6de6cd0847b;p=foosball

diff --git a/foorank.cpp b/foorank.cpp
index 096c185..f8ecf63 100644
--- a/foorank.cpp
+++ b/foorank.cpp
@@ -5,8 +5,9 @@
 #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;
@@ -74,19 +75,34 @@ double prodai(double a)
 // 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 /= 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);
-		double z = (x - mu2)/sigma2;
-		double gaussian = exp(-(z*z/2.0));
-		sum += step_size * probscore * gaussian;
+
+	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
@@ -221,24 +237,34 @@ void solve3x3(double *A, double *x, double *B)
 // of statistical moments.
 void estimate_musigma(vector<pair<double, double> > &curve, double &mu_result, double &sigma_result)
 {
-	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 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;
-		ex += (x1-x0) * xm * ym;
-		ex2 += (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;
 
-	ex /= sum_area;
-	ex2 /= sum_area;
+	area = (h/3.0) * area;
+	ex = (h/3.0) * ex / area;
+	ex2 = (h/3.0) * ex2 / area;
 
 	mu_result = ex;
 	sigma_result = sqrt(ex2 - ex * ex);
@@ -328,13 +354,13 @@ int main(int argc, char **argv)
 	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)));