3 * Simple linear regression
5 * (c) 2023 Straylight/Edgeware
8 /*----- Licensing notice --------------------------------------------------*
10 * This file is part of the mLib utilities library.
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.
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.
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,
28 /*----- Header files ------------------------------------------------------*/
35 /*----- Main code ---------------------------------------------------------*/
37 /* --- @linreg_init@ --- *
39 * Arguments: @struct linreg *lr@ = linear regression state
43 * Use: Initializes a linear-regression state ready for use.
46 void linreg_init(struct linreg
*lr
)
48 lr
->sum_x
= lr
->sum_x2
= lr
->sum_y
= lr
->sum_y2
= lr
->sum_x_y
= 0.0;
52 /* --- @linreg_update@ --- *
54 * Arguments: @struct linreg *lr@ = linear regression state
55 * @double x, y@ = point coordinates
59 * Use: Informs the linear regression machinery of a point.
61 * Note that the state size is constant, and independent of the
65 void linreg_update(struct linreg
*lr
, double x
, double y
)
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
;
72 /* --- @linreg_fit@ --- *
74 * Arguments: @const struct linreg *lr@ = linear regression state
75 * @double *m_out, *c_out, *r_out@ = where to write outputs
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.
86 * At least one point must have been given.
89 void linreg_fit(const struct linreg
*lr
,
90 double *m_out
, double *c_out
, double *r_out
)
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
;
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
;
100 cov_X_Y
= E_X_Y
- E_X
*E_Y
; var_X
= E_X2
- E2_X
; var_Y
= E_Y2
- E2_Y
;
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
);
108 /*----- That's all, folks -------------------------------------------------*/