Python:Fitting

This document contains examples of polynomial fitting, general linear regression, and nonlinear regression. In each section, there will be example code that may come in useful for later courses. The example code is based on the existence of a file in the same directory called Cantilever.dat that contains two columns of data - the first is an amount of mass (in kg) placed at the end of a beam and the second is a displacement, measured in inches, at the end of the beam. For EGR 103, this file is:

0.000000        0.005211
0.113510002     0.158707
0.227279999     0.31399
0.340790009     0.474619
0.455809998     0.636769
0.569320007     0.77989
0.683630005     0.936634
0.797140015     0.999986

Common Command Reference

All links below to NumPy v1.15 manual at NumPy v1.15 Manual; these commands show up in just about all the examples:

• numpy.ndarray.copy
• numpy.ndarray.min
• numpy.ndarray.max
• numpy.linspace
• numpy.loadtxt
• numpy.mean

Polynomial Fitting

Polynomial fits are those where the dependent data is related to some set of integer powers of the independent variable. MATLAB's built-in polyfit command can determine the coefficients of a polynomial fit.

Specific Command References

All links below to NumPy v1.15 manual at NumPy v1.15 Manual

• numpy.polyfit
• numpy.polyval

Example Code

In the example code below, n determines the order of the fit. Not much else would ever need to change.

 1 # %% Import modules
 2 import numpy as np
 3 import matplotlib.pyplot as plt
 4 
 5 
 6 # %% Load and manipulate data
 7 # Load data from Cantilever.dat
 8 beam_data = np.loadtxt('Cantilever.dat')
 9 # Copy data from each column into new variables
10 mass = beam_data[:, 0].copy()
11 disp = beam_data[:, 1].copy()
12 # Convert mass to force
13 force = mass * 9.81
14 # Convert disp to meters
15 disp = (disp * 2.54) / 100
16 
17 # %% Rename and create model data
18 x = force
19 y = disp
20 xmodel = np.linspace(np.min(x), np.max(x), 100)
21 
22 # %% Perform calculations
23 n = 1
24 p = np.polyfit(force, disp, n)
25 print(p)
26 
27 
28 
29 
30 # %% Generate estimates and model
31 yhat = np.polyval(p, x)
32 ymodel = np.polyval(p, xmodel)
33 
34 # %% Calculate statistics
35 st = np.sum((y - np.mean(y))**2)
36 sr = np.sum((y - yhat)**2)
37 r2 = (st - sr) / st
38 print('st: {}\nsr: {}\nr2: {}'.format(st, sr, r2))
39 
40 # %% Generate and save plots
41 plt.figure(1)
42 plt.clf()
43 plt.plot(x, y, 'ko', label='Data')
44 plt.plot(x, yhat, 'ks', label='Estimates', mfc='none')
45 plt.plot(xmodel, ymodel, 'k-', label='Model')
46 plt.grid(1)
47 plt.legend()

General Linear Regression

General linear regression involves finding some set of coefficients for fits that can be written as:

\( \hat{y}(x)=\sum_{j=1}^{M}a_j\phi_j(x) \)

where the \(a_j\) are the coefficients of the fit and the \(\phi_j\) are the specific functions of the independent variable that make up the fit.

Specific Command References

All links below to NumPy v1.15 manual at NumPy v1.15 Manual

• numpy.block
• numpy.linalg.lstsq
• numpy.reshape

Example Code

In the example code below, there is an example of a general linear fits of one variable. It is solving the same fit as given above, just in different way. Specifically it uses linear algebra to find the coefficients that minimize the sum of the squares of the estimate residuals for a general linear fit.

 1 # %% Import modules
 2 import numpy as np
 3 import matplotlib.pyplot as plt
 4 
 5 
 6 # %% Load and manipulate data
 7 # Load data from Cantilever.dat
 8 beam_data = np.loadtxt('Cantilever.dat')
 9 # Copy data from each column into new variables
10 mass = beam_data[:, 0].copy()
11 disp = beam_data[:, 1].copy()
12 # Convert mass to force
13 force = mass * 9.81
14 # Convert disp to meters
15 disp = (disp * 2.54) / 100
16 
17 # %% Rename and create model data
18 xv = np.reshape(force, (-1, 1))
19 yv = np.reshape(disp, (-1, 1))
20 xmodel = np.linspace(np.min(xv), np.max(xv), 100)
21 
22 # %% Perform calculations
23 def yfun(xe, coefs):
24     return coefs[0] * xe + coefs[1]
25 
26 a_mat = np.block([[xv**1, xv**0]])
27 pvec = np.linalg.lstsq(a_mat, yv, rcond=None)[0]
28 print(pvec)
29 
30 # %% Generate estimates and model
31 yhat = yfun(xv, pvec)
32 ymodel = yfun(xmodel, pvec)
33 
34 # %% Calculate statistics
35 st = np.sum((yv - np.mean(yv))**2)
36 sr = np.sum((yv - yhat)**2)
37 r2 = (st - sr) / st
38 print('st: {}\nsr: {}\nr2: {}'.format(st, sr, r2))
39 
40 # %% Generate and save plots
41 plt.figure(1)
42 plt.clf()
43 plt.plot(xv, yv, 'ko', label='Data')
44 plt.plot(xv, yhat, 'ks', label='Estimates', mfc='none')
45 plt.plot(xmodel, ymodel, 'k-', label='Model')
46 plt.grid(1)
47 plt.legend()

Nonlinear Regression

Nonlinear regression is both more powerful and more sensitive than linear regression. For inherently nonlinear fits, it will also produce a better \(S_r\) value than linearization since the nonlinear regression process is minimizing the \(S_r\) of the actual data rather than that of the transformed values. The sensitivity comes into play as the optimization routine may find local minima versus global minima. A good starting guess will work wonders.

Specific Command References

The link below is to the SciPy v1.1.0 reference guide at SciPy

• scipy.optimize.curve_fit

Example Code

Note in the example code that the initial guess gives 1 for the slope and 2 for the intercept. While these numbers are quite far from the optimized values of 0.0034 for the slope and 0.00055 for the intercept, the optimization routine is still able to find the correct value.

 1 # %% Import modules
 2 import numpy as np
 3 import matplotlib.pyplot as plt
 4 import scipy.optimize as opt
 5 
 6 # %% Load and manipulate data
 7 # Load data from Cantilever.dat
 8 beam_data = np.loadtxt('Cantilever.dat')
 9 # Copy data from each column into new variables
10 mass = beam_data[:, 0].copy()
11 disp = beam_data[:, 1].copy()
12 # Convert mass to force
13 force = mass * 9.81
14 # Convert disp to meters
15 disp = (disp * 2.54) / 100
16 
17 # %% Rename and create model data
18 x = force
19 y = disp
20 xmodel = np.linspace(np.min(x), np.max(x), 100)
21 
22 # %% Perform calculations
23 def yfun(x, *coefs):
24     return coefs[0] * x + coefs[1]
25 
26 popt = opt.curve_fit(yfun, x, y, [1, 2])[0]
27 print(popt)
28 
29 
30 # %% Generate estimates and model
31 yhat = yfun(x, *popt)
32 ymodel = yfun(xmodel, *popt)
33 
34 # %% Calculate statistics
35 st = np.sum((y - np.mean(y))**2)
36 sr = np.sum((y - yhat)**2)
37 r2 = (st - sr) / st
38 print('st: {}\nsr: {}\nr2: {}'.format(st, sr, r2))
39 
40 # %% Generate and save plots
41 plt.figure(1)
42 plt.clf()
43 plt.plot(x, y, 'ko', label='Data')
44 plt.plot(x, yhat, 'ks', label='Estimates', mfc='none')
45 plt.plot(xmodel, ymodel, 'k-', label='Model')
46 plt.grid(1)
47 plt.legend()