Linear Curve Fitting in IronPython QuickStart Sample

Illustrates how to fit linear combinations of curves to data using the LinearCurveFitter class and other classes in the Numerics.NET.Curves namespace in IronPython.

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

Overview

This QuickStart sample demonstrates how to perform linear least squares curve fitting using the LinearCurveFitter class in Numerics.NET.

The sample shows two practical applications of linear curve fitting. First, it demonstrates polynomial fitting using data from the NIST Statistical Reference Datasets library for calibrating load cells. The sample fits a quadratic polynomial to deflection measurements, showing both unweighted and weighted least squares approaches.

The second example shows how to fit arbitrary linear combinations of functions. It uses experimental data for the temperature-dependent thermal conductivity of copper. The theoretical model expresses conductance as a combination of inverse temperature and temperature squared terms. The sample demonstrates how to:

Create custom basis functions for fitting
Transform data to achieve linearity
Use the GeneralFunctionBasis class for arbitrary function combinations
Extract and interpret fitted parameters and their uncertainties
Calculate residuals to assess fit quality

The code includes examples of working with vectors, setting up fitting problems, handling weights, and interpreting statistical results. It serves as a practical introduction to linear least squares methods in Numerics.NET.

The code

import numerics

from System import Array
from System import Func

# The curve fitting classes reside in the 
# Extreme.Mathematics.Curves namespace.
from Extreme.Mathematics import *
from Extreme.Mathematics.Curves import *

# Illustrates least squares curve fitting of polynomials and
# other linear functions using the LinearCurveFitter class in the 
# Extreme.Mathematics.Curves namespace of the Extreme
# Optimization Numerical Libraries for .NET.

# This QuickStart sample illustrates linear least squares
# curve fitting using polynomials and linear combinations
# of arbitrary functions.

# Linear least squares fits are calculated using the
# LinearCurveFitter class:
fitter = LinearCurveFitter()

# We use data from the National Institute for Standards 
# and Technology's Statistical Reference Datasets library 
# at http:#www.itl.nist.gov/div898/strd/.

# Note that, due to round-off error, the results here will not be exactly
# the same as the NIST results, which were calculated using 500 digits
# of precision!

# We use the 'Pontius' dataset, which contains measurement data
# from the calibration of load cells. The independent variable is the load.
# The dependent variable is the deflection.
deflectionData = Vector([ .11019, .21956, .32949, .43899, .54803, .65694, \
    .76562, .87487, .98292, 1.09146, 1.20001, 1.30822, 1.41599, 1.52399, \
    1.63194, 1.73947, 1.84646, 1.95392, 2.06128, 2.16844, .11052, .22018, \
    .32939, .43886, .54798, .65739, .76596, .87474, .98300, 1.09150, \
    1.20004, 1.30818, 1.41613, 1.52408, 1.63159, 1.73965, 1.84696, \
    1.95445, 2.06177, 2.16829 ])
loadData =Vector([ 150, 300, 450, 600, 750, 900, 1050, 1200, 1350, 1500, \
    1650, 1800, 1950, 2100, 2250, 2400, 2550, 2700, 2850, 3000, 150, 300, \
    450, 600, 750, 900, 1050, 1200, 1350, 1500, 1650, 1800, 1950, 2100, \
    2250, 2400, 2550, 2700, 2850, 3000 ])

# You must supply the curve whose parameters will be
# fit to the data. The curve must inherit from LinearCombination.
#
# Here, we use a quadratic polynomial:
fitter.Curve = Polynomial(2)

# The X values go into the XValues property:
fitter.XValues = loadData
# ...and Y values go into the YValues property:
fitter.YValues = deflectionData

# 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
# The standard deviations associated with each parameter
# are available through the GetStandardDeviations method.
s = fitter.GetStandardDeviations()

print "Calibration of load cells"
print "    deflection = c1 + c2*load + c3*load^2 "
print "Solution:"
print "c1: {0:20.10e} {1:20.10e}".format(solution[0], s[0])
print "c2: {0:20.10e} {1:20.10e}".format(solution[1], s[1])
print "c3: {0:20.10e} {1:20.10e}".format(solution[2], s[2])

print "Residual sum of squares:", fitter.Residuals.Norm()

# 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.OneOverYSquared
# Refit the curve:
fitter.Fit()
solution = fitter.BestFitParameters
s = fitter.GetStandardDeviations()

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

#
# Fitting combinations of arbitrary functions
#

# The following example estimates the two parameters, c1 and c2
# in the theoretical model for conductance:
#     k(T) = 1 / (c1 / T + c2 * T*T)

temperature = Vector([ 12.2900, 13.7500, 14.8200, 16.1200, 18.0400, 18.6700, \
    20.5200, 22.6800, 25.1500, 27.7200, 30.2400, 33.2100, 36.4800, 39.8600, 50.4000 ])
conductance = Vector([ 25.3500, 27.8800, 29.9300, 30.4200, 31.0000, 31.9600, \
    32.4700, 30.3300, 31.1400, 27.4600, 23.2900, 20.7200, 17.2400, 14.7100,  9.5000 ])

# First, we transform the dependent variable:
y = Vector.Reciprocal(conductance)

# y is a linear combination of basis functions 1/T and T*T.
# Create a function basis object:
basisFunctions = Array[Func[float,float]]([ lambda x: 1 / x, lambda x: x**2 ])
basis = GeneralFunctionBasis(basisFunctions)
			
# Create a LinearCombination curve using this function basis:
curve = LinearCombination(basis)

# Set the curve fitter properties:
fitter.Curve = curve
fitter.XValues = temperature
fitter.YValues = y
# Reset the weights
fitter.WeightFunction = None
fitter.WeightVector = None

# Now compute the solution:
fitter.Fit()
solution = fitter.BestFitParameters
s = fitter.GetStandardDeviations()

# Print the results
print "Conductance of copper: k(T) = 1 / (c1/T + c2*T^2)"
print "Solution:"
print "c1: {0:20.10e} {1:20.10e}".format(solution[0], s[0])
print "c2: {0:20.10e} {1:20.10e}".format(solution[1], s[1])

print "Residual sum of squares:", fitter.Residuals.Norm()