Polynomial Least-Squares Fitting

pfit is a simple command-line program to do polynomial least-squares fitting to a given set of data, with uncertainties.

pfit.f90 is a corresponding Fortran 90 module to do the same fitting.

Synopsis

Let \( f(x) = a_0 + a_1 x + \ldots + a_d x^d \) be a polynomial of degree \( d \ge 0 \). The goal is to find the parameters \( a_0, \ldots, a_d \) such that \( f(x) \) "best" approximates a given set of data points \( (x_i, y_i, \sigma_i) \), \( i=1,\ldots,m \), where \( y_i \) has uncertainty \( \sigma_i \).

pfit finds the parameters \( a_0, \ldots, a_d\) such that the quantity \[ \chi^2 = \sum_{i=1}^{m} \frac{(f(x_i) - y_i)^2}{\sigma_i^2} \] is minimized.

Method

The simplest method to find \( a_0, \ldots, a_d\) is to set the derivative of \( \chi^2 \) to zero and derive a set of linear equations (the Normal Equations) that can be solved by simple matrix inversion. However, this method can be numerically unstable. This code instead uses a QR decomposition to minimize \( \chi^2 \), as described in Ref. [1].

Notes

This code requires LAPACK. The location of the LAPACK library can be specified in the Makefile.

Download

References

  1. N. Higham, Accuracy and Stability of Numerical Algorithms, Society for Industrial and Applied Mathematics (2002).