Commit | Line | Data |
---|---|---|
b64eb60f MW |
1 | /* -*-c-*- |
2 | * | |
3 | * Simple linear regression | |
4 | * | |
5 | * (c) 2023 Straylight/Edgeware | |
6 | */ | |
7 | ||
8 | /*----- Licensing notice --------------------------------------------------* | |
9 | * | |
10 | * This file is part of the mLib utilities library. | |
11 | * | |
12 | * mLib is free software: you can redistribute it and/or modify it under | |
13 | * the terms of the GNU Library General Public License as published by | |
14 | * the Free Software Foundation; either version 2 of the License, or (at | |
15 | * your option) any later version. | |
16 | * | |
17 | * mLib is distributed in the hope that it will be useful, but WITHOUT | |
18 | * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or | |
19 | * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Library General Public | |
20 | * License for more details. | |
21 | * | |
22 | * You should have received a copy of the GNU Library General Public | |
23 | * License along with mLib. If not, write to the Free Software | |
24 | * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, | |
25 | * USA. | |
26 | */ | |
27 | ||
28 | /*----- Header files ------------------------------------------------------*/ | |
29 | ||
30 | #include <assert.h> | |
31 | #include <math.h> | |
32 | ||
33 | #include "linreg.h" | |
34 | ||
35 | /*----- Main code ---------------------------------------------------------*/ | |
36 | ||
37 | /* --- @linreg_init@ --- * | |
38 | * | |
39 | * Arguments: @struct linreg *lr@ = linear regression state | |
40 | * | |
41 | * Returns: --- | |
42 | * | |
43 | * Use: Initializes a linear-regression state ready for use. | |
44 | */ | |
45 | ||
46 | void linreg_init(struct linreg *lr) | |
47 | { | |
48 | lr->sum_x = lr->sum_x2 = lr->sum_y = lr->sum_y2 = lr->sum_x_y = 0.0; | |
49 | lr->n = 0; | |
50 | } | |
51 | ||
52 | /* --- @linreg_update@ --- * | |
53 | * | |
54 | * Arguments: @struct linreg *lr@ = linear regression state | |
55 | * @double x, y@ = point coordinates | |
56 | * | |
57 | * Returns: --- | |
58 | * | |
59 | * Use: Informs the linear regression machinery of a point. | |
60 | * | |
61 | * Note that the state size is constant, and independent of the | |
62 | * number of points. | |
63 | */ | |
64 | ||
65 | void linreg_update(struct linreg *lr, double x, double y) | |
66 | { | |
67 | lr->sum_x += x; lr->sum_x2 += x*x; lr->sum_x_y += x*y; | |
68 | lr->sum_y += y; lr->sum_y2 += y*y; | |
69 | lr->n++; | |
70 | } | |
71 | ||
72 | /* --- @linreg_fit@ --- * | |
73 | * | |
74 | * Arguments: @const struct linreg *lr@ = linear regression state | |
75 | * @double *m_out, *c_out, *r_out@ = where to write outputs | |
76 | * | |
77 | * Returns: --- | |
78 | * | |
79 | * Use: Compute the best-fit line through the previously-specified | |
80 | * points. The line has the equation %$y = m x + c$%; %$m$% and | |
81 | * %$c$% are written to @*m_out@ and @*c_out@ respectively, and | |
82 | * the correlation coefficient %$r$% is written to @*r_out@. | |
83 | * Any (or all, but that would be useless) of the output | |
84 | * pointers may be null to discard that result. | |
85 | * | |
86 | * At least one point must have been given. | |
87 | */ | |
88 | ||
89 | void linreg_fit(const struct linreg *lr, | |
90 | double *m_out, double *c_out, double *r_out) | |
91 | { | |
92 | double E_X, E_X2, E2_X, E_X_Y, E_Y, E_Y2, E2_Y, n; | |
93 | double cov_X_Y, var_X, var_Y, m; | |
94 | assert(lr->n); | |
95 | ||
96 | n = lr->n; E_X_Y = lr->sum_x_y/n; | |
97 | E_X = lr->sum_x/n; E_X2 = lr->sum_x2/n; E2_X = E_X*E_X; | |
98 | E_Y = lr->sum_y/n; E_Y2 = lr->sum_y2/n; E2_Y = E_Y*E_Y; | |
99 | ||
100 | cov_X_Y = E_X_Y - E_X*E_Y; var_X = E_X2 - E2_X; var_Y = E_Y2 - E2_Y; | |
101 | ||
102 | m = cov_X_Y/var_X; | |
103 | if (m_out) *m_out = m; | |
104 | if (c_out) *c_out = E_Y - m*E_X; | |
105 | if (r_out) *r_out = cov_X_Y/sqrt(var_X*var_Y); | |
106 | } | |
107 | ||
108 | /*----- That's all, folks -------------------------------------------------*/ |