Surface fitting and more

So far, we have been dealing with nonlinear functions of one variable. In this section, we will show how to fit nonlinear functions of any number of variables of arbitrary types.

When defining what nonlinear fitting means for arbitrary inputs, the key observation is that in the definition of a nonlinear curve,

f(x) = f(a₁, a₂, ..., a_n; x) ,

the variable x can be of any type. The actual fitting algorithm only does calculations in terms of the curve parameters and in terms of functions defined by the user.

The NonlinearCurveFitter<T> class performs a nonlinear least squares fit where the type of the input is specified by the generic type argument, T.

The nonlinear function is defined by a delegate that maps a T and a Vector<T> value to a real number. The first argument is the X value. The second argument is the vector of function parameters. Likewise, the partial derivatives with respect to the curve parameters must be supplied as a delegate that maps a T value and a vector to a real vector.

The data is supplied as Vector<T> objects. The XValues and YValues properties specify the X and Y values of the data points, respectively. Note that in this case, the XValues vector has the element type of the generic type parameter.

By default, the fit is unweighted. All weights are set equal to 1. Weights can be specified in one of two ways. The first way is to set the WeightVector property to a vector with as many components as there are data points. Each component of the vector specifies the weight of the corresponding data point.

If error estimates for the Y values are available, then the GetWeightVectorFromErrors of the WeightFunctions class can be used to construct a vector with the weights that correspond to the errors. Its only argument is a Vector<T> that contains the error for each observation.

The Optimizer property gives access to the Levenberg-Marquardt optimizer that is used to compute the least squares solution. You can set tolerances and control convergence as set out in the chapter on optimization.

The Fit method performs the actual least squares calculation, and adjusts the parameters of the Curve to correspond to the least squares fit. The BestFitParameters property returns a Vector<T> containing the parameters of the curve. The GetStandardDeviations returns a Vector<T> containing the standard deviations of each parameter as estimated in the least squares calculation.

The Optimizer property also gives access to details of the calculation. It is useful to inspect its Status property to make sure the algorithm ended normally. A value other than Converged indicates some kind of problem with the algorithm.

Example: Fitting a Plane

We will use a very simple example: fitting a plane through a set of points.

Vector<Point> xValues = Vector.Create(
    new Point(-1, 1), new Point(1, 1),
    new Point(1, -1), new Point(-1, -1));
Vector<double> yValues = Vector.Create(0.0, 1.0, 2.0, 3.0);
NonlinearCurveFitter<Point> fitter = new NonlinearCurveFitter<Point>();
fitter.XValues = xValues;
fitter.YValues = yValues;
fitter.Function = (x, p) => p[0] + p[1] * x.X + p[2] * x.Y;
fitter.InitialGuess = Vector.Create<double>(3);
fitter.Fit();

Dim xValues = Vector.Create(
        New Point(-1, 1), New Point(1, 1),
        New Point(1, -1), New Point(-1, -1))
Dim yValues = Vector.Create(0.0, 1.0, 2.0, 3.0)
Dim fitter = New NonlinearCurveFitter(Of Point)()
fitter.XValues = xValues
fitter.YValues = yValues
fitter.Function = Function(x, p) p(0) + p(1) * x.X + p(2) * x.Y
fitter.InitialGuess = Vector.Create(Of Double)(3)
fitter.Fit()

No code example is currently available or this language may not be supported.

let xValues = Vector.Create(
                Point(-1.0, 1.0), Point(1.0, 1.0),
                Point(1.0, -1.0), Point(-1.0, -1.0))
let yValues = Vector.Create(0.0, 1.0, 2.0, 3.0)
let fitter = new NonlinearCurveFitter<Point>()
fitter.XValues <- xValues
fitter.YValues <- yValues
fitter.Function <- fun x p -> p.[0] + p.[1] * x.X + p.[2] * x.Y
fitter.InitialGuess <- Vector.Create<double>(3)
fitter.Fit()
//

The data points are defined here to be Point structures. The function is a simple linear function of the parameters.