lsfn.c File Reference

Detailed Description

Computes nth order polynomial least squares fit.

This file contains the implementation of the lsfn function, which computes a polynomial least squares fit of a specified order to given data points. It supports both weighted and unweighted fitting.

Author: Michael Borland, 1986.

Definition in file lsfn.c.

#include "matlib.h"
#include "mdb.h"

Go to the source code of this file.

Functions
int	p_merror (char *message)
long	lsfn (double *xd, double *yd, double *sy, long nd, long nf, double *coef, double *s_coef, double *chi, double *diff)
	Computes nth order polynomial least squares fit.

Function Documentation

lsfn()

long lsfn	(	double *	xd,
		double *	yd,
		double *	sy,
		long	nd,
		long	nf,
		double *	coef,
		double *	s_coef,
		double *	chi,
		double *	diff )

Computes nth order polynomial least squares fit.

This function performs an nth order polynomial least squares fit to the provided data. It supports both weighted and unweighted fitting based on the standard deviations provided.

Parameters

xd	Array of x data points.
yd	Array of y data points.
sy	Array of standard deviations of y data points. If NULL or all elements are equal, an unweighted fit is performed.
nd	Number of data points.
nf	Order of the polynomial fit (degree of the polynomial).
coef	Array to store the computed coefficients of the polynomial.
s_coef	Array to store the standard deviations of the coefficients.
chi	Pointer to store the reduced chi-squared value. If NULL, chi-squared is not computed.
diff	Array to store the differences between fitted and actual y values. If NULL, differences are not stored.

Returns

Returns 1 on success, or 0 on error.

Definition at line 34 of file lsfn.c.

  {
  long i, j, nt, unweighted;
  double xp, *x_i, x0;
  MATRIX *X = NULL, *Y = NULL, *Yp = NULL, *C = NULL, *C_1 = NULL, *Xt = NULL, *A = NULL, *Ca = NULL, *XtC = NULL, *XtCX = NULL, *T = NULL, *Tt = NULL, *TC = NULL;
  long status = 0;
 
  nt = nf + 1;
  if (nd < nt) {
    printf("error: insufficient data for requested order of fit\n");
    printf("(%ld data points, %ld terms in fit)\n", nd, nt);
    exit(1);
  }
 
  unweighted = 1;
  if (sy)
    for (i = 1; i < nd; i++)
      if (sy[i] != sy[0]) {
        unweighted = 0;
        break;
      }
 
  /* allocate matrices */
  m_alloc(&X, nd, nt);
  m_alloc(&Y, nd, 1);
  m_alloc(&Yp, nd, 1);
  m_alloc(&Xt, nt, nd);
  if (!unweighted) {
    m_alloc(&C, nd, nd);
    m_alloc(&C_1, nd, nd);
    m_zero(C);
    m_zero(C_1);
  }
  m_alloc(&A, nt, 1);
  m_alloc(&Ca, nt, nt);
  m_alloc(&XtC, nt, nd);
  m_alloc(&XtCX, nt, nt);
  m_alloc(&T, nt, nd);
  m_alloc(&Tt, nd, nt);
  m_alloc(&TC, nt, nd);
 
  /* Compute X, Y, C, C_1.  X[i][j] = (xd[i])^j. Y[i][0] = yd[i].
   * C   = delta(i,j)*sy[i]^2  (covariance matrix of yd)
   * C_1 = INV(C)
   */
  for (i = 0; i < nd; i++) {
    x_i = X->a[i];
    x0 = xd[i];
    xp = 1.0;
    Y->a[i][0] = yd[i];
    if (!unweighted) {
      C->a[i][i] = sqr(sy[i]);
      C_1->a[i][i] = 1 / C->a[i][i];
    }
    for (j = 0; j < nt; j++) {
      x_i[j] = xp;
      xp *= x0;
    }
  }
 
  /* Compute A, the matrix of coefficients.
   * Weighted least-squares solution is A = INV(Xt.INV(C).X).Xt.INV(C).y 
   * Unweighted solution is A = INV(Xt.X).Xt.y 
   */
  if (unweighted) {
    /* eliminating 2 matrix operations makes this much faster than a weighted fit
       * if there are many data points.
       */
    if (!m_trans(Xt, X))
      { status = p_merror("transposing X"); goto cleanup; }
    if (!m_mult(XtCX, Xt, X))
      { status = p_merror("multiplying Xt.X"); goto cleanup; }
    if (!m_invert(XtCX, XtCX))
      { status = p_merror("inverting XtCX"); goto cleanup; }
    if (!m_mult(T, XtCX, Xt))
      { status = p_merror("multiplying XtX.Xt"); goto cleanup; }
    if (!m_mult(A, T, Y))
      { status = p_merror("multiplying T.Y"); goto cleanup; }
 
    /* Compute covariance matrix of A, Ca = (T.Tt)*C[0][0] */
    if (!m_trans(Tt, T))
      { status = p_merror("computing transpose of T"); goto cleanup; }
    if (!m_mult(Ca, T, Tt))
      { status = p_merror("multiplying T.Tt"); goto cleanup; }
    if (!m_scmul(Ca, Ca, sy ? sqr(sy[0]) : 1))
      { status = p_merror("multiplying T.Tt by scalar"); goto cleanup; }
  } else {
    if (!m_trans(Xt, X))
      { status = p_merror("transposing X"); goto cleanup; }
    if (!m_mult(XtC, Xt, C_1))
      { status = p_merror("multiplying Xt.C_1"); goto cleanup; }
    if (!m_mult(XtCX, XtC, X))
      { status = p_merror("multiplying XtC.X"); goto cleanup; }
    if (!m_invert(XtCX, XtCX))
      { status = p_merror("inverting XtCX"); goto cleanup; }
    if (!m_mult(T, XtCX, XtC))
      { status = p_merror("multiplying XtCX.XtC"); goto cleanup; }
    if (!m_mult(A, T, Y))
      { status = p_merror("multiplying T.Y"); goto cleanup; }
 
    /* Compute covariance matrix of A, Ca = T.C.Tt */
    if (!m_mult(TC, T, C))
      { status = p_merror("multiplying T.C"); goto cleanup; }
    if (!m_trans(Tt, T))
      { status = p_merror("computing transpose of T"); goto cleanup; }
    if (!m_mult(Ca, TC, Tt))
      { status = p_merror("multiplying TC.Tt"); goto cleanup; }
  }
 
  for (i = 0; i < nt; i++) {
    coef[i] = A->a[i][0];
    if (s_coef)
      s_coef[i] = sqrt(Ca->a[i][i]);
  }
 
  /* Compute Yp = X.A, use to compute chi-squared */
  if (chi) {
    if (!m_mult(Yp, X, A))
      { status = p_merror("multiplying X.A"); goto cleanup; }
    *chi = 0;
    for (i = 0; i < nd; i++) {
      xp = (Yp->a[i][0] - yd[i]);
      if (diff != NULL)
        diff[i] = xp;
      xp /= sy ? sy[i] : 1;
      *chi += xp * xp;
    }
    if (nd != nt)
      *chi /= (nd - nt);
  }
 
  status = 1;
 
cleanup:
  m_free(&X);
  m_free(&Y);
  m_free(&Yp);
  m_free(&Xt);
  if (!unweighted) {
    m_free(&C);
    m_free(&C_1);
  }
  m_free(&A);
  m_free(&Ca);
  m_free(&XtC);
  m_free(&XtCX);
  m_free(&T);
  m_free(&Tt);
  m_free(&TC);
  return (status);
}