Nonlinear Curve Fitting in F# QuickStart Sample

Illustrates nonlinear least squares curve fitting of predefined and user-defined curves using the NonlinearCurveFitter class in F#.

View this sample in: C# Visual Basic IronPython

// Illustrates nonlinear least squares curve fitting using the
// NonlinearCurveFitter class in the 
// Numerics.NET.Curves namespace of Numerics.NET.

#light

open System

open Numerics.NET
// The curve fitting classes reside in the 
// Numerics.NET.Curves namespace.
open Numerics.NET.Curves
// The predefined non-linear curves reside in the 
// Numerics.NET.Curves namespace.
open Numerics.NET.Curves.Nonlinear
// Vectors reside in the Numerics.NET.Mathemaics.LinearAlgebra
// namespace
open Numerics.NET.LinearAlgebra
// The non-linear least squares optimizer resides in the
// Numerics.NET.Optimization namespace.
open Numerics.NET.Optimization

// The license is verified at runtime. We're using a 30 day trial key here.
// For more information, see:
//     https://numerics.net/trial-key
let licensed = Numerics.NET.License.Verify("64542-18980-57619-62268")

// Nonlinear least squares fits are calculated using the
// NonlinearCurveFitter class:
let fitter = NonlinearCurveFitter()

// In the first example, we fit a dose response curve
// to a data set that includes error information.

// The data points must be supplied as Vector objects:
let dose = 
    Vector.Create(1.46247, 2.3352, 4.0, 
        7.0, 12.0, 18.0, 23.0, 30.0, 40.0, 60.0, 90.0, 160.0, 290.0, 490.0, 860.0)
let response = 
    Vector.Create(95.49073, 95.14551, 94.86448, 
        92.66762, 85.36377, 74.72183, 62.76747, 51.04137, 38.20257, 
        28.01712, 19.40086, 13.18117, 9.87161, 7.64622, 7.21826)
let error = 
    Vector.Create(4.74322, 4.74322, 4.74322, 
        4.63338, 4.26819, 3.73609, 3.13837, 3.55207, 3.91013, 
        2.40086, 2.6, 3.65906, 2.49358, 2.38231, 2.36091)

// You must supply the curve whose parameters will be
// fit to the data. The curve must inherit from NonlinearCurve.
// The FourParameterLogistic curve is one of several
// predefined nonlinear curves:
let doseResponseCurve = FourParameterLogisticCurve()

// Now we set the curve fitter's Curve property:
fitter.Curve <- doseResponseCurve
// The GetInitialFitParameters method sets the curve parameters
// to initial values appropriate for the data:
fitter.InitialGuess <- doseResponseCurve.GetInitialFitParameters(dose, response)

// and the data values:
fitter.XValues <- dose
fitter.YValues <- response
// The GetWeightVectorFromErrors method of the WeightFunctions
// class lets us convert the error values to weights:
fitter.WeightVector <- WeightFunctions.GetWeightVectorFromErrors(error)

// The Fit method performs the actual calculation.
let result = fitter.Fit()
// The standard deviations associated with each parameter
// are available through the GetStandardDeviations method.
let s = fitter.GetStandardDeviations()

// We can now print the results:
printfn "Dose response curve"

printfn "Initial value: %10.6f +/- %.4f" doseResponseCurve.InitialValue s.[0]
printfn "Final value:   %10.6f +/- %.4f" doseResponseCurve.FinalValue s.[1]
printfn "Center:        %10.6f +/- %.4f" doseResponseCurve.Center s.[2]
printfn "Hill slope:    %10.6f +/- %.4f" doseResponseCurve.HillSlope s.[3]

// We can also show some statistics about the calculation:
printfn "Residual sum of squares: %A" (fitter.Residuals.Norm())
// The Optimizer property returns the MultidimensionalOptimization object
// used to perform the calculation:
printfn "# iterations: %d" fitter.Optimizer.IterationsNeeded
printfn "# function evaluations: %d" fitter.Optimizer.EvaluationsNeeded

printfn ""

//
// Defining your own nonlinear curve
//

// In this example, we use one of the datasets (MGH10) 
// from the National Institute for Statistics and Technology
// (NIST) Statistical Reference Datasets.
// See http://www.itl.nist.gov/div898/strd for details

let fitter2 = NonlinearCurveFitter()
// Here, we need to define our own curve.
// The MyCurve class is defined below.

// This is our nonlinear curve implementation. For details, see 
// http://www.itl.nist.gov/div898/strd/nls/data/mgh10.shtml
// You must inherit from NonlinearCurve:
type MyCurve() as this =
    // Call the base constructor with the number of parameters.
    inherit NonlinearCurve(3)

    do
        // It is convenient to set common starting values
        // for the curve parameters in the constructor:
        this.Parameters.[0] <- 0.2
        this.Parameters.[1] <- 40000.0
        this.Parameters.[2] <- 2500.0

    override this.ValueAt x =
        this.Parameters.[0] * exp(this.Parameters.[1] / (x + this.Parameters.[2]))
        
    override this.SlopeAt x =
        this.Parameters.[0] * this.Parameters.[1] * exp(this.Parameters.[1] / (x + this.Parameters.[2]))
            / (pown (x + this.Parameters.[2]) 2)

    // The FillPartialDerivatives evaluates the partial derivatives
    // with respect to the curve parameters, and returns
    // the result in a vector. If you don't supply this method, 
    // a numerical approximation is used.
    override this.FillPartialDerivatives (x, f) =
        let exp = Math.Exp(this.Parameters.[1] / (x + this.Parameters.[2]))
        f.[0] <- exp
        f.[1] <- this.Parameters.[0] * exp / (x + this.Parameters.[2])
        f.[2] <- -this.Parameters.[0] * this.Parameters.[1] * exp / (pown (x + this.Parameters.[2]) 2)

fitter2.Curve <- MyCurve()

// The data is provided as Vector objects.
// X values go into the XValues property...
fitter2.XValues <- Vector.Create(
    [|
        5.000000E+01; 5.500000E+01; 6.000000E+01; 6.500000E+01;
        7.000000E+01; 7.500000E+01; 8.000000E+01; 8.500000E+01;
        9.000000E+01; 9.500000E+01; 1.000000E+02; 1.050000E+02;
        1.100000E+02; 1.150000E+02; 1.200000E+02; 1.250000E+02    
    |])
// ...and Y values go into the YValues property:
fitter2.YValues <- Vector.Create(
    [|
        3.478000E+04; 2.861000E+04; 2.365000E+04; 1.963000E+04;
        1.637000E+04; 1.372000E+04; 1.154000E+04; 9.744000E+03;
        8.261000E+03; 7.030000E+03; 6.005000E+03; 5.147000E+03;
        4.427000E+03; 3.820000E+03; 3.307000E+03; 2.872000E+03
    |])
fitter2.InitialGuess <- Vector.Create(fitter2.Curve.Parameters.ToArray())
fitter2.WeightVector <- null
// The Fit method performs the actual calculation:
let result2 = fitter2.Fit()

// A Vector containing the parameters of the best fit
// can be obtained through the
// BestFitParameters property.
let solution2 = fitter2.BestFitParameters
let s2 = fitter2.GetStandardDeviations()

printfn "NIST Reference Data Set"
printfn "Solution:"
printfn "b1: %20f %20f" solution2.[0] s2.[0]
printfn "b2: %20f %20f" solution2.[1] s2.[1]
printfn "b3: %20f %20f" solution2.[2] s2.[2]

printfn "Certified values:"
printfn "b1: %20f %20f" 5.6096364710E-03 1.5687892471E-04 
printfn "b2: %20f %20f" 6.1813463463E+03 2.3309021107E+01 
printfn "b3: %20f %20f" 3.4522363462E+02 7.8486103508E-01 

// Now let's redo the same operation, but with observations weighted
// by 1/Y^2. To do this, we set the WeightFunction property.
// The WeightFunctions class defines a set of ready-to-use weight functions.
fitter2.WeightFunction <- WeightFunctions.OneOverX
// Refit the curve:
let result3 = fitter2.Fit()
let solution3 = fitter2.BestFitParameters
let s3 = fitter2.GetStandardDeviations()

// The solution is slightly different:
printfn "Solution (weighted observations):"
printfn "b1: %20f %20f" solution3.[0] s3.[0]
printfn "b2: %20f %20f" solution3.[1] s3.[1]
printfn "b3: %20f %20f" solution3.[2] s3.[2]

printf "Press Enter key to exit..."
Console.ReadLine() |> ignore