@@@ wip

[mLib] / utils / linreg.c

/* -*-c-*-
 *
 * Simple linear regression
 *
 * (c) 2023 Straylight/Edgeware
 */

/*----- Licensing notice --------------------------------------------------*
 *
 * This file is part of the mLib utilities library.
 *
 * mLib is free software: you can redistribute it and/or modify it under
 * the terms of the GNU Library General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or (at
 * your option) any later version.
 *
 * mLib is distributed in the hope that it will be useful, but WITHOUT
 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
 * FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Library General Public
 * License for more details.
 *
 * You should have received a copy of the GNU Library General Public
 * License along with mLib.  If not, write to the Free Software
 * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307,
 * USA.
 */

/*----- Header files ------------------------------------------------------*/

#include <assert.h>
#include <math.h>

#include "linreg.h"

/*----- Main code ---------------------------------------------------------*/

/* --- @linreg_init@ --- *
 *
 * Arguments:	@struct linreg *lr@ = linear regression state
 *
 * Returns:	---
 *
 * Use:		Initializes a linear-regression state ready for use.
 */

void linreg_init(struct linreg *lr)
{
  lr->sum_x = lr->sum_x2 = lr->sum_y = lr->sum_y2 = lr->sum_x_y = 0.0;
  lr->n = 0;
}

/* --- @linreg_update@ --- *
 *
 * Arguments:	@struct linreg *lr@ = linear regression state
 *		@double x, y@ = point coordinates
 *
 * Returns:	---
 *
 * Use:		Informs the linear regression machinery of a point.
 *
 *		Note that the state size is constant, and independent of the
 *		number of points.
 */

void linreg_update(struct linreg *lr, double x, double y)
{
  lr->sum_x += x; lr->sum_x2 += x*x; lr->sum_x_y += x*y;
  lr->sum_y += y; lr->sum_y2 += y*y;
  lr->n++;
}

/* --- @linreg_fit@ --- *
 *
 * Arguments:	@const struct linreg *lr@ = linear regression state
 *		@double *m_out, *c_out, *r_out@ = where to write outputs
 *
 * Returns:	---
 *
 * Use:		Compute the best-fit line through the previously-specified
 *		points.  The line has the equation %$y = m x + c$%; %$m$% and
 *		%$c$% are written to @*m_out@ and @*c_out@ respectively, and
 *		the correlation coefficient %$r$% is written to @*r_out@.
 *		Any (or all, but that would be useless) of the output
 *		pointers may be null to discard that result.
 *
 *		At least one point must have been given.
 */

void linreg_fit(const struct linreg *lr,
		double *m_out, double *c_out, double *r_out)
{
  double E_X, E_X2, E2_X, E_X_Y, E_Y, E_Y2, E2_Y, n;
  double cov_X_Y, var_X, var_Y, m;
  assert(lr->n);

  n = lr->n; E_X_Y = lr->sum_x_y/n;
  E_X = lr->sum_x/n; E_X2 = lr->sum_x2/n; E2_X = E_X*E_X;
  E_Y = lr->sum_y/n; E_Y2 = lr->sum_y2/n; E2_Y = E_Y*E_Y;

  cov_X_Y = E_X_Y - E_X*E_Y; var_X = E_X2 - E2_X; var_Y = E_Y2 - E2_Y;

  m = cov_X_Y/var_X;
  if (m_out) *m_out = m;
  if (c_out) *c_out = E_Y - m*E_X;
  if (r_out) *r_out = cov_X_Y/sqrt(var_X*var_Y);
}

/*----- That's all, folks -------------------------------------------------*/