Polynomial Curve Best Fit

Home → Knowledgebase → Delphi → Polynomial Curve Best Fit

Background

We recently had a requirement for a customer to calculate a best-fit curve for a set of points. The software we were implementing is used to drive a metalworking machine and it had to produce an output which was a smooth curve, based on a set of sampled points. Having worked on similar problems before, I had a rough idea how this should work mathematically but as is usually the case when developing software, I was sure that someone else would have already solved this particular problem.

Requirements

We have a set of points with X and Y coordinates. The software needs to calculate a curve which is a best fit for the points. In doing so, it should also introduce additional intermediate points so that the curve is as smooth as possible. For example, for an input set of 400 points, the software will output 5,000 points which all lie on the curve of best fit.

Research

A quick Google search identified some possible Delphi solutions to this problem.

The one that looked most promising was the least squares curve fit unit by David Taylor: https://www.satsignal.eu/software/components.html#CurveFit

This was a single Delphi unit that had no dependencies on third-party components and simply used the System.Math unit.
The routine was written in Delphi 5, but the test project compiled in the latest version of Delphi without any issues, as one would expect with a purely mathematical library.

Using the Library

Inputs

The main function accepts:

• An array of floating-point X values.
• An array of floating-point Y values.
• The number of terms to use.

The number of terms relates to the degree of the polynomial formula which will be used.

A single term equates to a horizontal line on a graph, e.g.   y = 2

The use of two terms equates to a straight line on a graph, e.g.   y = 3x + 2

Using three terms equates to a quadratic equation, e.g.   y = 4x² + 3x + 2

Outputs

The calculation function will populate an array of floating-point values which correspond to the coefficients of the equation. For example, if we want to use three terms (a quadratic equation), we can pass an array of three floating-point numbers. In the quadratic example above, the coefficients would be 2, 3 and 4 respectively. The calculation also returns a correlation coefficient, which tells us how closely the curve correlates with the input points provided. The closer this coefficient is to 1, the more accurate the fit.

Results

When using the library, unit tests were written using a set of test points provided by our customer. The points provided were typical of the samples they would be using in production. With these sample points, the most accurate curve was produced when calling the function with six terms. This equates to a polynomial function of the fifth degree.

When using more than six terms, there were occasional floating-point overflow errors produced, presumably because the coefficients calculated where too small to be represented in a Double type in Delphi. When using fewer than six terms, the routine still calculated a curve of best fit, but of course the fit was not as accurate. With fewer terms, the correlation coefficient was lower.

With six terms, we received six coefficients back in the array provided to the function.

Once we have the coefficients, we can calculate the Y value for any value of X. To achieve the desired curve, an array of 5,000 evenly spaced X points was calculated, bounded by the maximum and minimum X values of the sample.

function CalculateY(const X: Double; const Coefficients: TArray<Double>): Double;
begin
  var yc := 0.0;
  var xc := 1.0;

  for var i := Low(Coefficients) to High(Coefficients) do
  begin
    yc := yc + Coefficients[i] * xc;
    xc := xc * X;
  end;

  Result := yc;
end;

Here is a sample of some of the input points:

In the software, the resulting best fit curve is calculated and overlaid over the points as follows:

Implementation

The function uses a least-squares curve fitting algorithm using Gauss-Jordan elimination.

This is a well known algorithm for finding a best fit polynomial curve for a set of points.

Conclusion

We have seen that a mathematical library originally written in Turbo Pascal and early versions of Delphi can be used in the latest versions of Delphi without modification.

This library can produce accurate and smooth curves for the best fit of a set of data points.
The higher the number of terms and hence the degree of the polynomial, the more accurate the fit.