Nonlinear Curve Fitting in IronPython QuickStart Sample

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

This sample is also available in: C#, Visual Basic, F#.

Overview

This QuickStart sample demonstrates how to perform nonlinear curve fitting using the NonlinearCurveFitter class in Numerics.NET.

The sample shows two main approaches to nonlinear curve fitting:

  1. Using a predefined curve type (FourParameterLogisticCurve) to fit a dose-response curve to data that includes measurement errors. This demonstrates how to:
    • Set up the curve fitter with x/y data points and measurement errors
    • Convert errors to weights for weighted least squares fitting
    • Obtain fit parameters and their standard deviations
    • Access diagnostic information about the fitting process
  2. Creating and fitting a custom nonlinear curve using both manual implementation and automatic differentiation:
    • Shows how to create a custom curve class by inheriting from NonlinearCurve
    • Demonstrates implementing value, slope and partial derivative calculations
    • Shows how to use automatic differentiation as an alternative to manual implementation
    • Uses the NIST Statistical Reference Dataset MGH10 to validate the implementation
    • Demonstrates both unweighted and weighted fitting using built-in weight functions

The sample includes extensive error handling and provides examples of accessing diagnostic information such as residuals, iteration counts, and function evaluations.

The code

import numerics

from math import exp

# The curve fitting classes reside in the 
# Extreme.Mathematics.Curves namespace.
from Extreme.Mathematics.Curves import *
# The predefined non-linear curves reside in the 
# Extreme.Mathematics.Curves namespace.
from Extreme.Mathematics.Curves.Nonlinear import *
# Vectors reside in the Extreme.Mathemaics namespace
from Extreme.Mathematics import *
# The non-linear least squares optimizer resides in the
# Extreme.Mathematics.Optimization namespace.
from Extreme.Mathematics.Optimization import *

# Illustrates nonlinear least squares curve fitting using the
# NonlinearCurveFitter class in the 
# Extreme.Mathematics.Curves namespace of the Extreme
# Optimization Numerical Libraries for .NET.

# Nonlinear least squares fits are calculated using the
# NonlinearCurveFitter class:
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:
dose = Vector.Create[float](1.46247, 2.3352, 4, 7, 12, 18, 23, 30, \
    40, 60, 90, 160, 290, 490, 860)
response = Vector.Create[float](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)
error = Vector.Create[float](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:
doseResponseCurve = FourParameterLogisticCurve()


# Now we set the curve fitter's Curve property:
fitter.Curve = doseResponseCurve
# and the data values:
fitter.XValues = dose
fitter.YValues = response
# The GetInitialFitParameters method returns a set of
# initial values appropriate for the data:
fitter.InitialGuess = doseResponseCurve.GetInitialFitParameters(dose, 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.
fitter.Fit()
# The standard deviations associated with each parameter
# are available through the GetStandardDeviations method.
s = fitter.GetStandardDeviations()

# We can now print the results:
print "Dose response curve"

print "Initial value: {0:10.6f} +/- {1:.4f}".format(doseResponseCurve.InitialValue, s[0])
print "Final value:   {0:10.6f} +/- {1:.4f}".format(doseResponseCurve.FinalValue, s[1])
print "Center:        {0:10.6f} +/- {1:.4f}".format(doseResponseCurve.Center, s[2])
print "Hill slope:    {0:10.6f} +/- {1:.4f}".format(doseResponseCurve.HillSlope, s[3])

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

print 

#
# 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

fitter = NonlinearCurveFitter()

from symbolic import *
(x,a,b,c) = symbols(['x','a','b','c'])

fitter.Curve = NonlinearCurve.FromExpression(fun([x,a,b,c],  a * exp(b / (x + c))))

# In this case, we have to specify the initial values
# for the parameters:
fitter.InitialGuess = Vector.Create[float](0.2, 40000, 2500)

# The data is provided as Vector objects.
# X values go into the XValues property...
fitter.XValues = Vector.Create[float]( \
	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:
fitter.YValues = Vector.Create[float]( \
	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)
fitter.WeightVector = None
# The Fit method performs the actual calculation:
fitter.Fit()

# A Vector containing the parameters of the best fit
# can be obtained through the
# BestFitParameters property.
solution = fitter.BestFitParameters
s = fitter.GetStandardDeviations()

print "NIST Reference Data Set"
print "Solution:"
print "b1: {0:20.10e} {1:20.10e}".format(solution[0], s[0])
print "b2: {0:20.10e} {1:20.10e}".format(solution[1], s[1])
print "b3: {0:20.10e} {1:20.10e}".format(solution[2], s[2])

print "Certified values:"
print "b1: {0:20.10e} {1:20.10e}".format(5.6096364710E-03, 1.5687892471E-04)
print "b2: {0:20.10e} {1:20.10e}".format(6.1813463463E+03, 2.3309021107E+01)
print "b3: {0:20.10e} {1:20.10e}".format(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.
fitter.WeightFunction = WeightFunctions.OneOverX
# Refit the curve:
fitter.Fit()
solution = fitter.BestFitParameters
s = fitter.GetStandardDeviations()

# The solution is slightly different:
print "Solution (weighted observations):"
print "b1: {0:20.10e} {1:20.10e}".format(solution[0], s[0])
print "b2: {0:20.10e} {1:20.10e}".format(solution[1], s[1])
print "b3: {0:20.10e} {1:20.10e}".format(solution[2], s[2])