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.

View this sample in: C# Visual Basic F#

```Python
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])
```