chiark - git - mdw - mLib/blob - utils/linreg.c

   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 -------------------------------------------------*/