-
Notifications
You must be signed in to change notification settings - Fork 0
/
linreg.h
103 lines (83 loc) · 2.62 KB
/
linreg.h
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
#pragma once
struct LinearFunction {
double m;
double b;
/* correlation coeff */
double r;
constexpr double operator()(double x) const { return m * x + b; }
};
namespace LinearRegression {
struct Point {
double x;
double y;
};
constexpr double square(double x) { return x * x; }
constexpr double sqrtNewtonRaphson(double x, double curr, double prev) {
return curr == prev ? curr
: sqrtNewtonRaphson(x, 0.5 * (curr + x / curr), curr);
}
constexpr double sqrt(double x) {
return x < 0 || x >= INFINITY ? NAN : sqrtNewtonRaphson(x, x, 0);
}
class Sums {
constexpr double sumx(int n, Point const points[]) {
return n == 0 ? 0 : points[n - 1].x + sumx(n - 1, points);
}
constexpr double sumx2(int n, Point const points[]) {
return n == 0 ? 0 : square(points[n - 1].x) + sumx2(n - 1, points);
}
constexpr double sumy(int n, Point const points[]) {
return n == 0 ? 0 : points[n - 1].y + sumy(n - 1, points);
}
constexpr double sumy2(int n, Point const points[]) {
return n == 0 ? 0 : square(points[n - 1].y) + sumy2(n - 1, points);
}
constexpr double sumxy(int n, Point const points[]) {
return n == 0 ? 0
: points[n - 1].x * points[n - 1].y + sumxy(n - 1, points);
}
public:
double x;
double x2;
double y;
double y2;
double xy;
constexpr Sums(int n, Point const points[])
: x(sumx(n, points)),
x2(sumx2(n, points)),
y(sumy(n, points)),
y2(sumy2(n, points)),
xy(sumxy(n, points)) {}
};
struct Result {
bool ok;
LinearFunction f;
static constexpr Result success(double m, double b, double r) {
return {true, {m, b, r}};
}
static constexpr Result fail() { return {false}; }
};
constexpr Result compute1(int n, Sums const &sums, double denom) {
return Result::success(
/* m */ (n * sums.xy - sums.x * sums.y) / denom,
/* b */ (sums.y * sums.x2 - sums.x * sums.xy) / denom,
/* r */ (sums.xy - sums.x * sums.y / n) /
sqrt((sums.x2 - square(sums.x) / n) *
(sums.y2 - square(sums.y) / n)));
}
constexpr Result decide(int n, Sums const &sums, double denom) {
// If denominator is 0, we have a singular matrix => can't solve the problem.
return denom == 0 ? Result::fail() : compute1(n, sums, denom);
}
constexpr Result compute(int n, Sums const &sums) {
return decide(n, sums, n * sums.x2 - square(sums.x));
}
class Line {
Result result_;
public:
template <int N>
constexpr Line(Point const (&points)[N])
: result_(compute(N, Sums(N, points))) {}
constexpr operator LinearFunction() const { return result_.f; }
};
} // namespace LinearRegression