@@@ man wip

[mLib] / utils / linreg.h

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

#ifndef MLIB_LINREG_H
#define MLIB_LINREG_H

#ifdef __cplusplus
  extern "C" {
#endif

/* Theory.  (This should be well-known.)
 *
 * We have a number of data points (x_i, y_i) for 0 <= i < n.  We believe
 * they lie close to a straight line y = m x + c, for unknown constants m and
 * c.  Let e_i = m x_i + c - y_i.  We seek the parameters which minimize
 * E = ∑_{0\le i<n} e_i^2.
 *
 * We begin by simplifying
 *
 *	E = ∑ e_i^2
 *	  = ∑ (m x_i + c - y_i)^2
 *	  = ∑ (m^2 x_i^2 + c^2 + y_i^2 + 2 c m x_i - 2 m x_i y_i - 2 c y_i) .
 *
 * (where all sums are over 0 <= i < n).  Taking partial derivatives,
 *
 *	∂E/∂m = ∑ (2 m x_i^2 + 2 c x_i - 2 x_i y_i)
 *	      = 2 (m ∑ x_i^2 + c ∑ x_i - ∑ x_i y_i)
 *
 * and
 *
 *	∂E/∂c = ∑ (2 c + 2 m x_i - 2 y_i)
 *	      = 2 (n c + m ∑ x_i - ∑ y_i) .
 *
 * Now we solve for m and c such that the partial derivatives vanish.  Note
 * that the second partial derivatives are nonnegative, being 2 ∑ x_i^2 and 2
 * n respectively, so we shall indeed find a minimum.
 *
 * Restating the problem,
 *
 *	m ∑ x_i^2 + c ∑ x_i - ∑ x_i y_i = 0
 *	m ∑ x_i   + c n     - ∑ y_i     = 0 .
 *
 * Writing E[X] = (∑ x_i)/n, E[X Y] = (∑ x_i y_i)/n, etc., this gives
 *
 *	m n E[X^2] + c n E[X] - n E[X Y] = 0
 *	m n E[X]   + c n      - n E[Y]   = 0 .
 *
 * We see a common factor of n, so we can take that out.  Now multiply the
 * latter equation by E[X] and subtract, giving
 *
 *	m (E[X^2] - E[X]^2) - (E[X Y] - E[X] E[Y])
 *
 * whence
 *
 *	m = (E[X Y] - E[X] E[Y])/(E[X^2] - E[X]^2) .
 *
 * The numerator is known as the covariance of X and Y, written
 *
 *	cov(X, Y) = σ_{XY} = E[X Y] - E[X] E[Y] ;
 *
 * the denominator is the variance of X, written
 *
 *	var(X) = σ_X = E[X^2] - E[X]^2 .
 *
 * The second equation gives
 *
 *	c = E[Y] - m E[X] .
 *
 * Also of interest is the correlation coefficient
 *
 *	r = σ_{XY}/√(σ_X σ_Y) ,
 *
 * which I'm not going to derive here.
 */

/*----- Data structures ---------------------------------------------------*/

struct linreg {
  double sum_x, sum_x2, sum_y, sum_y2, sum_x_y;
  unsigned long n;
};
#define LINREG_INIT { 0.0, 0.0, 0.0, 0.0, 0.0, 0 }

/*----- Functions provided ------------------------------------------------*/

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

extern void linreg_init(struct linreg */*lr*/);

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

extern void linreg_update(struct linreg */*lr*/, double /*x*/, double /*y*/);

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

extern void linreg_fit(const struct linreg */*lr*/,
		       double */*m_out*/, double */*c_out*/,
		       double */*r_out*/);

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

#ifdef __cplusplus
  }
#endif

#endif