Difference between revisions of "Python:Fitting"

From PrattWiki
Jump to navigation Jump to search
 
(33 intermediate revisions by the same user not shown)
Line 1: Line 1:
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 <code>Cantilever.dat</code> 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:  
+
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 <code>Cantilever.dat</code> 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:
 
<source lang="matlab">
 
<source lang="matlab">
 
0.000000        0.005211
 
0.000000        0.005211
Line 23: Line 23:
 
== Common Code ==
 
== Common Code ==
 
The processes for loading and manipulating the original data, calculating the statistics, and creating the graph are the same across the examples.  To simplify codes, those processes should be contained in a separate file called <code>fitting_common.py</code> as given below:
 
The processes for loading and manipulating the original data, calculating the statistics, and creating the graph are the same across the examples.  To simplify codes, those processes should be contained in a separate file called <code>fitting_common.py</code> as given below:
<source lang=python>
+
<syntaxhighlight lang=python line>
 
import numpy as np
 
import numpy as np
 
import matplotlib.pyplot as plt
 
import matplotlib.pyplot as plt
  
def get_beam_data():  
+
 
 +
def get_beam_data():
 
     # Load data from Cantilever.dat
 
     # Load data from Cantilever.dat
     beam_data = np.loadtxt('Cantilever.dat')
+
     beam_data = np.loadtxt("Cantilever.dat")
 
     # Copy data from each column into new variables
 
     # Copy data from each column into new variables
 
     mass = beam_data[:, 0].copy()
 
     mass = beam_data[:, 0].copy()
Line 37: Line 38:
 
     # Convert disp to meters
 
     # Convert disp to meters
 
     disp = (disp * 2.54) / 100
 
     disp = (disp * 2.54) / 100
   
+
 
 
     fmodel = np.linspace(np.min(force), np.max(force), 100)
 
     fmodel = np.linspace(np.min(force), np.max(force), 100)
   
+
 
 
     # Return values
 
     # Return values
 
     return force, disp, fmodel
 
     return force, disp, fmodel
 +
  
 
def calc_stats(y, yhat, to_print=1):
 
def calc_stats(y, yhat, to_print=1):
 
     # Calculate statistics
 
     # Calculate statistics
     st = np.sum((y - np.mean(y))**2)
+
     st = np.sum((y - np.mean(y)) ** 2)
     sr = np.sum((y - yhat)**2)
+
     sr = np.sum((y - yhat) ** 2)
 
     r2 = (st - sr) / st
 
     r2 = (st - sr) / st
 
     if to_print:
 
     if to_print:
         print('st: {:.3e}\nsr: {:.3e}\nr2: {:.3e}'.format(st, sr, r2))
+
         print("st: {:.3e}\nsr: {:.3e}\nr2: {:.3e}".format(st, sr, r2))
 
     return st, sr, r2
 
     return st, sr, r2
 +
  
 
def make_plot(x, y, yhat, xmodel, ymodel, fig_num=1):
 
def make_plot(x, y, yhat, xmodel, ymodel, fig_num=1):
    plt.figure(fig_num).clf()
+
     fig, ax = plt.subplots(num=fig_num, clear=True)
     fig, ax = plt.subplots(num=fig_num)
+
     ax.plot(x, y, "ko", label="Data")
     ax.plot(x, y, 'ko', label='Data')
+
     ax.plot(x, yhat, "ks", label="Estimates", mfc="none")
     ax.plot(x, yhat, 'ks', label='Estimates', mfc='none')
+
     ax.plot(xmodel, ymodel, "k-", label="Model")
     ax.plot(xmodel, ymodel, 'k-', label='Model')
 
 
     ax.grid(True)
 
     ax.grid(True)
     ax.legend(location='best')
+
     ax.legend(loc="best")
 
     fig.tight_layout()
 
     fig.tight_layout()
</source>
+
    return fig, ax
 +
 
 +
</syntaxhighlight>
  
 
== Polynomial Fitting ==
 
== Polynomial Fitting ==
Line 73: Line 77:
 
=== Example Code ===
 
=== Example Code ===
 
In the example code below, <code>n</code> determines the order of the fit.  Not much else would ever need to change.
 
In the example code below, <code>n</code> determines the order of the fit.  Not much else would ever need to change.
<syntaxhighlight  lang='python' line highlight='8-20'>
+
<syntaxhighlight  lang='python' line highlight='9-23'>
 
# %% Import modules
 
# %% Import modules
 
import numpy as np
 
import numpy as np
 
from fitting_common import *
 
from fitting_common import *
 +
  
 
# %% Load and manipulate data
 
# %% Load and manipulate data
Line 85: Line 90:
 
p = np.polyfit(x, y, n)
 
p = np.polyfit(x, y, n)
 
print(p)
 
print(p)
 +
 +
  
  
Line 105: Line 112:
 
General linear regression involves finding some set of coefficients for fits that can be written as:
 
General linear regression involves finding some set of coefficients for fits that can be written as:
 
<center><math>
 
<center><math>
\hat{y}(x)=\sum_{j=1}^{M}c_j\phi_j(x)
+
\hat{y}(x)=\sum_{m=0}^{M-1}a_m\phi_m(x)
 
</math></center>
 
</math></center>
where the <math>c_j</math> are the coefficients of the fit and the <math>\phi_j</math> are the specific functions of the independent variable that make up the fit.
+
where the <math>a_m</math> are the coefficients of the fit and the <math>\phi_m(x)</math> are the specific functions of the independent variable that make up the fit.
 
=== Specific Command References ===
 
=== Specific Command References ===
 
All links below to NumPy v1.15 manual at [https://docs.scipy.org/doc/numpy-1.15.0/index.html NumPy v1.15 Manual]
 
All links below to NumPy v1.15 manual at [https://docs.scipy.org/doc/numpy-1.15.0/index.html NumPy v1.15 Manual]
Line 117: Line 124:
 
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.  In this code, variables ending in "v" explicitly need to be column vectors while variables ending in "e" can either be 1D arrays or 2D arrays (or lists).
 
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.  In this code, variables ending in "v" explicitly need to be column vectors while variables ending in "e" can either be 1D arrays or 2D arrays (or lists).
  
Note - if you get an error about "TypeError: must be real number, not NoneType" then on line 27, replace None with -1.
+
Note - if you get an error about "TypeError: must be real number, not NoneType" then on line 15, replace None with -1.
<syntaxhighlight  lang='python' line highlight='8-22>
+
<syntaxhighlight  lang='python' line highlight=9-23>
 
# %% Import modules
 
# %% Import modules
 
import numpy as np
 
import numpy as np
 
from fitting_common import *
 
from fitting_common import *
 +
  
 
# %% Load and manipulate data
 
# %% Load and manipulate data
Line 129: Line 137:
 
def yfun(xe, coefs):
 
def yfun(xe, coefs):
 
     return coefs[0] * xe**1 + coefs[1] * xe**0
 
     return coefs[0] * xe**1 + coefs[1] * xe**0
 +
 +
 
# Reshape data for block matrices
 
# Reshape data for block matrices
 
xv = np.reshape(x, (-1, 1))
 
xv = np.reshape(x, (-1, 1))
 
yv = np.reshape(y, (-1, 1))
 
yv = np.reshape(y, (-1, 1))
a_mat = np.block([[xv**1, xv**0]])
+
phi_mat = np.block([[xv**1, xv**0]])
pvec = np.linalg.lstsq(a_mat, yv, rcond=None)[0]
+
pvec = np.linalg.lstsq(phi_mat, yv, rcond=None)[0]
 
print(pvec)
 
print(pvec)
  
Line 145: Line 155:
 
# %% Generate plots
 
# %% Generate plots
 
make_plot(x, y, yhat, xmodel, ymodel)
 
make_plot(x, y, yhat, xmodel, ymodel)
 +
</syntaxhighlight>
 +
 +
=== Multidimensional General Linear Regression ===
 +
The processes above can be used for data sets with more than one independent variable.  The key is that anything involving linear algebra needs to be done with column vector versions of your data while anything graphical needs to be done with matrix versions.  The statistics can be done with columns if using <code>sum(thing)</code> or with either columns or matrices if using <code>thing.sum()</code>.    The sum command used on a 2-D array adds up the column; the sum method applied to an array adds up all the entries.
 +
 +
Multidimensional general linear regression involves finding some set of coefficients for fits that can be written as:
 +
<center><math>
 +
\hat{y}(x_1, x_2, x_3, ...)=\sum_{m=0}^{M-1}a_m\phi_m(x_1, x_2, x_3, ...)
 +
</math></center>
 +
where the <math>a_m</math> are the coefficients of the fit and the <math>\phi_m(x)</math> are the specific functions of the independent variables $$x_1, x_2, x_3, ...$$ that make up the fit.  If specifically looking at a function of two independent variables, this may be written as:
 +
<center><math>
 +
\hat{z}(x, y)=\sum_{m=0}^{M-1}a_m\phi_m(x, y)
 +
</math></center>
 +
where now $$x$$ and $$y$$ are independent variables and $$z$$ is the dependent variable.
 +
 +
==== Example Code ====
 +
In the example code below, there is an example of a general linear fits of two variables.  Specifically, the code will try to find the best general second order model to fit a collection of random numbers.  In this code, variables ending in "v" explicitly need to be column vectors, variables ending in "m" explicitly need to be matrices, and variables ending in "e" can either be 1D arrays or 2D arrays (or lists).
 +
 +
<syntaxhighlight  lang='python'>
 +
import numpy as np
 +
from fitting_common import *
 +
import matplotlib.pyplot as plt
 +
from mpl_toolkits.mplot3d import axes3d
 +
from matplotlib import cm
 +
 +
# %% Create data
 +
xm, ym = np.meshgrid(np.arange(-5, 5.1, 5), np.arange(-5, 5.1, 5))
 +
zm = np.random.uniform(size=xm.shape)
 +
xmodelm, ymodelm = np.meshgrid(np.linspace(xm.min(), xm.max(), 21), np.linspace(ym.min(), ym.max(), 21))
 +
 +
# %% Perform calculations
 +
def zfun(xe, ye, coefs):
 +
    return coefs[0] * xe**0 + coefs[1] * xe + coefs[2] * ye + coefs[3] * xe**2 + coefs[4] * xe * ye + coefs[5] * ye**2
 +
 +
 +
# Reshape data for block matrices
 +
xv = np.reshape(xm, (-1, 1))
 +
yv = np.reshape(ym, (-1, 1))
 +
zv = np.reshape(zm, (-1, 1))
 +
phi_mat = np.block([[xv**0, xv, yv, xv**2, xv*yv, yv**2]])
 +
pvec = np.linalg.lstsq(phi_mat, zv, rcond=None)[0]
 +
print(pvec)
 +
 +
# %% Generate estimates and model
 +
zhatv = zfun(xv, yv, pvec) # for stats
 +
zhatm = zfun(xm, ym, pvec) # for graphics
 +
zmodelm = zfun(xmodelm, ymodelm, pvec)
 +
 +
# %% Calculate statistics
 +
calc_stats(zv, zhatv, 1)
 +
 +
# %% Generate plots
 +
# Make model plot
 +
fig = plt.figure(num=1, clear=True)
 +
ax = fig.add_subplot(1, 1, 1, projection='3d')
 +
 +
surfhandle = ax.plot_surface(xmodelm, ymodelm, zmodelm, cmap=cm.autumn)
 +
ax.scatter(xm, ym, zm)
 +
ax.set(xlabel='x', ylabel='y', zlabel='z', title='Second Order Model')
 +
fig.colorbar(surfhandle)
 +
fig.tight_layout()
 +
 +
# Make error plot
 +
fig = plt.figure(num=2, clear=True)
 +
ax = fig.add_subplot(1, 1, 1, projection='3d')
 +
 +
ax.plot_surface(xm, ym, zm-zhatm, cmap=cm.autumn)
 +
ax.set(xlabel='x', ylabel='y', zlabel='$z-\hat{z}$', title='Second Order Model Error')
 +
fig.tight_layout()
 +
 
</syntaxhighlight>
 
</syntaxhighlight>
  
Line 151: Line 231:
  
 
=== Specific Command References ===
 
=== Specific Command References ===
The link below is to the SciPy v1.1.0 reference guide at [https://docs.scipy.org/doc/scipy/reference/index.html SciPy]
+
The link below is to the SciPy reference guide at [https://docs.scipy.org/doc/scipy/reference/index.html SciPy]
 
* [https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.curve_fit.html scipy.optimize.curve_fit]
 
* [https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.curve_fit.html scipy.optimize.curve_fit]
  
=== Example Code ===
+
=== Steps for Nonlinear Regression ===
Note in the example code that the initial guess gives 0.6 for the slope and 0.1 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.  That is not always the case - try to find an initial guess close to the actual answer.
+
The critical parts of solving for the nonlinear regression involve defining the function, setting the initial conditions, and understanding the output from the opt.curve_fit command.
<syntaxhighlight  lang='python' line highlight='8-20'>
+
 
 +
==== Defining The Function ====
 +
When the opt.curve_fit command calls the function to find the best values for the coefficients, it calls it with an independent argument for each of the coefficients.  That is to say, if the command is looking to solve a function for a straight line which has two coefficients - a slope and an intercept - the opt.curve_fit command will call the function with ''three'' inputs - the independent values, a guess for the slope, and a guess for the intercept.  This is different from the paradigm we have used so far with a single list containing all the coefficients.  In order to make the code more modular, we will use *args to take the coefficients in the function call and store them all in a single list.  For example, to define a function for a straight line model using non-linear regression, the function call would be:
 +
<syntaxhighlight  lang='python'>
 +
def yfun(x, *coefs):
 +
    return coefs[0] * x**1 + coefs[1] * x**0
 +
</syntaxhighlight>
 +
Now the function can be called two different ways: either with three inputs or with two inputs where the first is the value or set of values for the independent variable and the second is an unpacked list or tuple with the items in it.  That is to say, the following will calculate the <math>y</math> values for a line with the <math>x</math> values being a linearly spaced array of 5 values between 0 and 2 given a slope of 10 and an intercept of 8:
 +
<syntaxhighlight  lang='python'>
 +
In [1]: yfun(np.linspace(0, 2, 5), 10, 8)
 +
Out[1]: array([ 8., 13., 18., 23., 28.])
 +
 
 +
In [2]: yfun(np.linspace(0, 2, 5), *[10, 8])
 +
Out[2]: array([ 8., 13., 18., 23., 28.])
 +
</syntaxhighlight>
 +
The * before the list is key to unpack it into multiple parts; those multiple parts are then packed into the *arg.
 +
 
 +
==== Setting the Initial Conditions ====
 +
Unlike the general linear regression method, which will find the best coefficients for a linear fit without needing any initial guess, the nonlinear regression method requires a good initial guess.  Bad initial guesses can lead to errors or incorrect / non-optimal solutions.  Note in the example code below that the initial guess gives 0.6 for the slope and 0.1 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.  That is not always the case - try to find an initial guess close to the actual answer. Sometimes the coefficients have names that may help - if some constant is called <code>a_min</code> you should perhaps guess the minimum <code>a</code> value in your data set.
 +
 
 +
==== Understanding the Output ====
 +
The opt.curve_fit command returns two items in a tuple: the parameters themselves and some statistical information.  Since you only want the first of these, it makes sense to put a [0] at the end of the command to just grab the parameter values. Remember that you will still need to unpack the list of parameters when you call your function.
 +
 
 +
=== Example for a Straight Line ===
 +
 
 +
<syntaxhighlight  lang='python' line highlight=9-23>
 
# %% Import modules
 
# %% Import modules
 +
import numpy as np
 +
from fitting_common import *
 
import scipy.optimize as opt
 
import scipy.optimize as opt
from fitting_common import *
 
  
 
# %% Load and manipulate data
 
# %% Load and manipulate data
Line 166: Line 272:
 
# %% Perform calculations
 
# %% Perform calculations
 
def yfun(x, *coefs):
 
def yfun(x, *coefs):
     return coefs[0] * x**1 + coefs[1] * x**0
+
     return coefs[0] * x ** 1 + coefs[1] * x ** 0
 +
 
  
 
popt = opt.curve_fit(yfun, x, y, [0.6, 0.1])[0]
 
popt = opt.curve_fit(yfun, x, y, [0.6, 0.1])[0]
 
print(popt)
 
print(popt)
 +
  
  
Line 195: Line 303:
 
</math></center>
 
</math></center>
 
the code is
 
the code is
<syntaxhighlight lang=python line start=8 highlight="3,7">
+
<syntaxhighlight lang=python line start=9 highlight="3,7">
 
# %% Perform calculations
 
# %% Perform calculations
 
def yfun(xe, coefs):
 
def yfun(xe, coefs):
     return coefs[0] * xe**1 + coefs[1] * xe**0
+
     return coefs[0] * xe ** 1 + coefs[1] * xe ** 0
 
# Reshape data for block matrices
 
# Reshape data for block matrices
 
xv = np.reshape(x, (-1, 1))
 
xv = np.reshape(x, (-1, 1))
 
yv = np.reshape(y, (-1, 1))
 
yv = np.reshape(y, (-1, 1))
a_mat = np.block([[xv**1, xv**0]])
+
phi_mat = np.block([[xv ** 1, xv ** 0]])
pvec = np.linalg.lstsq(a_mat, yv, rcond=None)[0]
+
pvec = np.linalg.lstsq(phi_mat, yv, rcond=None)[0]
 
print(pvec)
 
print(pvec)
 
</syntaxhighlight>
 
</syntaxhighlight>
Line 215: Line 323:
 
</math></center>
 
</math></center>
 
then the code would change to:
 
then the code would change to:
<syntaxhighlight lang=python line start=8 highlight="3,7">
+
<syntaxhighlight lang=python line start=9 highlight="3,7">
 
# %% Perform calculations
 
# %% Perform calculations
 
def yfun(xe, coefs):
 
def yfun(xe, coefs):
     return coefs[0] * np.cos(xe) + coefs[1] * np.sin(xe) + coefs[2] * xe**0
+
     return coefs[0] * np.cos(xe) + coefs[1] * np.sin(xe) + coefs[2] * xe ** 0
 
# Reshape data for block matrices
 
# Reshape data for block matrices
 
xv = np.reshape(x, (-1, 1))
 
xv = np.reshape(x, (-1, 1))
 
yv = np.reshape(y, (-1, 1))
 
yv = np.reshape(y, (-1, 1))
a_mat = np.block([[np.cos(xv), np.sin(xv), xv**0]])
+
phi_mat = np.block([[np.cos(xv), np.sin(xv), xv ** 0]])
pvec = np.linalg.lstsq(a_mat, yv, rcond=None)[0]
+
pvec = np.linalg.lstsq(phi_mat, yv, rcond=None)[0]
 
print(pvec)
 
print(pvec)
 
</syntaxhighlight>
 
</syntaxhighlight>
Line 234: Line 342:
 
</math></center>
 
</math></center>
 
the code is
 
the code is
<syntaxhighlight lang=python line start=8 highlight="3,5">
+
<syntaxhighlight lang=python line start=9 highlight="3,5">
 
# %% Perform calculations
 
# %% Perform calculations
 
def yfun(x, *coefs):
 
def yfun(x, *coefs):
     return coefs[0] * x**1 + coefs[1] * x**0
+
     return coefs[0] * x ** 1 + coefs[1] * x ** 0
  
 
popt = opt.curve_fit(yfun, x, y, [0.6, 0.1])[0]
 
popt = opt.curve_fit(yfun, x, y, [0.6, 0.1])[0]
Line 251: Line 359:
 
</math></center>
 
</math></center>
 
then the code would change to:
 
then the code would change to:
<syntaxhighlight lang=python line start=8 highlight="3,5">
+
<syntaxhighlight lang=python line start=9 highlight="3,5">
 
# %% Perform calculations
 
# %% Perform calculations
 
def yfun(x, *coefs):
 
def yfun(x, *coefs):
     return coefs[0] * np.cos(coefs[1]*x + coefs[2]) + coefs[3] * x**0
+
     return coefs[0] * np.cos(coefs[1]*x + coefs[2]) + coefs[3] * x ** 0
  
 
popt = opt.curve_fit(yfun, x, y, [.01, 1, np.pi/4, .0125])[0]
 
popt = opt.curve_fit(yfun, x, y, [.01, 1, np.pi/4, .0125])[0]

Latest revision as of 02:33, 21 November 2022

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:

In the scripts below, common code has a regular background while code that differs from script to script will be highlighted in yellow.

Common Code

The processes for loading and manipulating the original data, calculating the statistics, and creating the graph are the same across the examples. To simplify codes, those processes should be contained in a separate file called fitting_common.py as given below:

 1 import numpy as np
 2 import matplotlib.pyplot as plt
 3 
 4 
 5 def get_beam_data():
 6     # Load data from Cantilever.dat
 7     beam_data = np.loadtxt("Cantilever.dat")
 8     # Copy data from each column into new variables
 9     mass = beam_data[:, 0].copy()
10     disp = beam_data[:, 1].copy()
11     # Convert mass to force
12     force = mass * 9.81
13     # Convert disp to meters
14     disp = (disp * 2.54) / 100
15 
16     fmodel = np.linspace(np.min(force), np.max(force), 100)
17 
18     # Return values
19     return force, disp, fmodel
20 
21 
22 def calc_stats(y, yhat, to_print=1):
23     # Calculate statistics
24     st = np.sum((y - np.mean(y)) ** 2)
25     sr = np.sum((y - yhat) ** 2)
26     r2 = (st - sr) / st
27     if to_print:
28         print("st: {:.3e}\nsr: {:.3e}\nr2: {:.3e}".format(st, sr, r2))
29     return st, sr, r2
30 
31 
32 def make_plot(x, y, yhat, xmodel, ymodel, fig_num=1):
33     fig, ax = plt.subplots(num=fig_num, clear=True)
34     ax.plot(x, y, "ko", label="Data")
35     ax.plot(x, yhat, "ks", label="Estimates", mfc="none")
36     ax.plot(xmodel, ymodel, "k-", label="Model")
37     ax.grid(True)
38     ax.legend(loc="best")
39     fig.tight_layout()
40     return fig, ax

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

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 from fitting_common import *
 4 
 5 
 6 # %% Load and manipulate data
 7 x, y, xmodel = get_beam_data()
 8 
 9 # %% Perform calculations
10 n = 1
11 p = np.polyfit(x, y, n)
12 print(p)
13 
14 
15 
16 
17 
18 
19 
20 
21 # %% Generate estimates and model
22 yhat = np.polyval(p, x)
23 ymodel = np.polyval(p, xmodel)
24 
25 # %% Calculate statistics
26 calc_stats(y, yhat, 1)
27 
28 # %% Generate plots
29 make_plot(x, y, yhat, xmodel, ymodel)

General Linear Regression

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

\( \hat{y}(x)=\sum_{m=0}^{M-1}a_m\phi_m(x) \)

where the \(a_m\) are the coefficients of the fit and the \(\phi_m(x)\) 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

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. In this code, variables ending in "v" explicitly need to be column vectors while variables ending in "e" can either be 1D arrays or 2D arrays (or lists).

Note - if you get an error about "TypeError: must be real number, not NoneType" then on line 15, replace None with -1.

 1 # %% Import modules
 2 import numpy as np
 3 from fitting_common import *
 4 
 5 
 6 # %% Load and manipulate data
 7 x, y, xmodel = get_beam_data()
 8 
 9 # %% Perform calculations
10 def yfun(xe, coefs):
11     return coefs[0] * xe**1 + coefs[1] * xe**0
12 
13 
14 # Reshape data for block matrices
15 xv = np.reshape(x, (-1, 1))
16 yv = np.reshape(y, (-1, 1))
17 phi_mat = np.block([[xv**1, xv**0]])
18 pvec = np.linalg.lstsq(phi_mat, yv, rcond=None)[0]
19 print(pvec)
20 
21 # %% Generate estimates and model
22 yhat = yfun(x, pvec)
23 ymodel = yfun(xmodel, pvec)
24 
25 # %% Calculate statistics
26 calc_stats(y, yhat, 1)
27 
28 # %% Generate plots
29 make_plot(x, y, yhat, xmodel, ymodel)

Multidimensional General Linear Regression

The processes above can be used for data sets with more than one independent variable. The key is that anything involving linear algebra needs to be done with column vector versions of your data while anything graphical needs to be done with matrix versions. The statistics can be done with columns if using sum(thing) or with either columns or matrices if using thing.sum(). The sum command used on a 2-D array adds up the column; the sum method applied to an array adds up all the entries.

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

\( \hat{y}(x_1, x_2, x_3, ...)=\sum_{m=0}^{M-1}a_m\phi_m(x_1, x_2, x_3, ...) \)

where the \(a_m\) are the coefficients of the fit and the \(\phi_m(x)\) are the specific functions of the independent variables $$x_1, x_2, x_3, ...$$ that make up the fit. If specifically looking at a function of two independent variables, this may be written as:

\( \hat{z}(x, y)=\sum_{m=0}^{M-1}a_m\phi_m(x, y) \)

where now $$x$$ and $$y$$ are independent variables and $$z$$ is the dependent variable.

Example Code

In the example code below, there is an example of a general linear fits of two variables. Specifically, the code will try to find the best general second order model to fit a collection of random numbers. In this code, variables ending in "v" explicitly need to be column vectors, variables ending in "m" explicitly need to be matrices, and variables ending in "e" can either be 1D arrays or 2D arrays (or lists).

import numpy as np
from fitting_common import *
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import axes3d
from matplotlib import cm

# %% Create data
xm, ym = np.meshgrid(np.arange(-5, 5.1, 5), np.arange(-5, 5.1, 5))
zm = np.random.uniform(size=xm.shape)
xmodelm, ymodelm = np.meshgrid(np.linspace(xm.min(), xm.max(), 21), np.linspace(ym.min(), ym.max(), 21))

# %% Perform calculations
def zfun(xe, ye, coefs):
    return coefs[0] * xe**0 + coefs[1] * xe + coefs[2] * ye + coefs[3] * xe**2 + coefs[4] * xe * ye + coefs[5] * ye**2


# Reshape data for block matrices
xv = np.reshape(xm, (-1, 1))
yv = np.reshape(ym, (-1, 1))
zv = np.reshape(zm, (-1, 1))
phi_mat = np.block([[xv**0, xv, yv, xv**2, xv*yv, yv**2]])
pvec = np.linalg.lstsq(phi_mat, zv, rcond=None)[0]
print(pvec)

# %% Generate estimates and model
zhatv = zfun(xv, yv, pvec) # for stats
zhatm = zfun(xm, ym, pvec) # for graphics
zmodelm = zfun(xmodelm, ymodelm, pvec)

# %% Calculate statistics
calc_stats(zv, zhatv, 1)

# %% Generate plots
# Make model plot
fig = plt.figure(num=1, clear=True)
ax = fig.add_subplot(1, 1, 1, projection='3d')

surfhandle = ax.plot_surface(xmodelm, ymodelm, zmodelm, cmap=cm.autumn)
ax.scatter(xm, ym, zm)
ax.set(xlabel='x', ylabel='y', zlabel='z', title='Second Order Model')
fig.colorbar(surfhandle)
fig.tight_layout()

# Make error plot
fig = plt.figure(num=2, clear=True)
ax = fig.add_subplot(1, 1, 1, projection='3d')

ax.plot_surface(xm, ym, zm-zhatm, cmap=cm.autumn)
ax.set(xlabel='x', ylabel='y', zlabel='$z-\hat{z}$', title='Second Order Model Error')
fig.tight_layout()

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 reference guide at SciPy

Steps for Nonlinear Regression

The critical parts of solving for the nonlinear regression involve defining the function, setting the initial conditions, and understanding the output from the opt.curve_fit command.

Defining The Function

When the opt.curve_fit command calls the function to find the best values for the coefficients, it calls it with an independent argument for each of the coefficients. That is to say, if the command is looking to solve a function for a straight line which has two coefficients - a slope and an intercept - the opt.curve_fit command will call the function with three inputs - the independent values, a guess for the slope, and a guess for the intercept. This is different from the paradigm we have used so far with a single list containing all the coefficients. In order to make the code more modular, we will use *args to take the coefficients in the function call and store them all in a single list. For example, to define a function for a straight line model using non-linear regression, the function call would be:

def yfun(x, *coefs):
    return coefs[0] * x**1 + coefs[1] * x**0

Now the function can be called two different ways: either with three inputs or with two inputs where the first is the value or set of values for the independent variable and the second is an unpacked list or tuple with the items in it. That is to say, the following will calculate the \(y\) values for a line with the \(x\) values being a linearly spaced array of 5 values between 0 and 2 given a slope of 10 and an intercept of 8:

In [1]: yfun(np.linspace(0, 2, 5), 10, 8)
Out[1]: array([ 8., 13., 18., 23., 28.])

In [2]: yfun(np.linspace(0, 2, 5), *[10, 8])
Out[2]: array([ 8., 13., 18., 23., 28.])

The * before the list is key to unpack it into multiple parts; those multiple parts are then packed into the *arg.

Setting the Initial Conditions

Unlike the general linear regression method, which will find the best coefficients for a linear fit without needing any initial guess, the nonlinear regression method requires a good initial guess. Bad initial guesses can lead to errors or incorrect / non-optimal solutions. Note in the example code below that the initial guess gives 0.6 for the slope and 0.1 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. That is not always the case - try to find an initial guess close to the actual answer. Sometimes the coefficients have names that may help - if some constant is called a_min you should perhaps guess the minimum a value in your data set.

Understanding the Output

The opt.curve_fit command returns two items in a tuple: the parameters themselves and some statistical information. Since you only want the first of these, it makes sense to put a [0] at the end of the command to just grab the parameter values. Remember that you will still need to unpack the list of parameters when you call your function.

Example for a Straight Line

 1 # %% Import modules
 2 import numpy as np
 3 from fitting_common import *
 4 import scipy.optimize as opt
 5 
 6 # %% Load and manipulate data
 7 x, y, xmodel = get_beam_data()
 8 
 9 # %% Perform calculations
10 def yfun(x, *coefs):
11     return coefs[0] * x ** 1 + coefs[1] * x ** 0
12 
13 
14 popt = opt.curve_fit(yfun, x, y, [0.6, 0.1])[0]
15 print(popt)
16 
17 
18 
19 
20 
21 # %% Generate estimates and model
22 yhat = yfun(x, *popt)
23 ymodel = yfun(xmodel, *popt)
24 
25 # %% Calculate statistics
26 calc_stats(y, yhat, 1)
27 
28 # %% Generate plots
29 make_plot(x, y, yhat, xmodel, ymodel)

Example changes for different models

Polynomial

For the polynomial fitting model, really the only thing that would change would be the order of the fit and thus the value of \(n\) on line 23 of that code.

General Linear

For the general linear fit, the two places things will change will be in the function definition on line 10 and in the creation of the block matrix on line 14; the changes will mirror each other. For the straight line model - formally:

\( \hat{y}(x)=mx^1+bx^0 \)

the code is

 9 # %% Perform calculations
10 def yfun(xe, coefs):
11     return coefs[0] * xe ** 1 + coefs[1] * xe ** 0
12 # Reshape data for block matrices
13 xv = np.reshape(x, (-1, 1))
14 yv = np.reshape(y, (-1, 1))
15 phi_mat = np.block([[xv ** 1, xv ** 0]])
16 pvec = np.linalg.lstsq(phi_mat, yv, rcond=None)[0]
17 print(pvec)

where pvec[0][0] will be the slope \(m\) and pvec[1][0] will be the intercept \(b\). If the model were changed to, say,

\( \hat{y}(x)=a\cos(x)+b\sin(x)+c \)

which, formally, is

\( \hat{y}(x)=a\cos(x)+b\sin(x)+cx^0 \)

then the code would change to:

 9 # %% Perform calculations
10 def yfun(xe, coefs):
11     return coefs[0] * np.cos(xe) + coefs[1] * np.sin(xe) + coefs[2] * xe ** 0
12 # Reshape data for block matrices
13 xv = np.reshape(x, (-1, 1))
14 yv = np.reshape(y, (-1, 1))
15 phi_mat = np.block([[np.cos(xv), np.sin(xv), xv ** 0]])
16 pvec = np.linalg.lstsq(phi_mat, yv, rcond=None)[0]
17 print(pvec)

and pvec will be a 2D array with three rows.

Nonlinear

For the nonlinear fit, the two places things will change will be in the function definition on line 10 and in the initial parameter value guess on line 12; there need to be as many values in the initial parameter list as there are parameters. For the straight line model - formally:

\( \hat{y}(x)=mx^1+bx^0 \)

the code is

 9 # %% Perform calculations
10 def yfun(x, *coefs):
11     return coefs[0] * x ** 1 + coefs[1] * x ** 0
12 
13 popt = opt.curve_fit(yfun, x, y, [0.6, 0.1])[0]
14 print(popt)

where popt[0] will be the slope \(m\) and popt[1] will be the intercept \(b\). If the model were changed to, say,

\( \hat{y}(x)=a\cos(bx+c)+d \)

which, formally, is

\( \hat{y}(x)=a\cos(bx+c)+dx^0 \)

then the code would change to:

 9 # %% Perform calculations
10 def yfun(x, *coefs):
11     return coefs[0] * np.cos(coefs[1]*x + coefs[2]) + coefs[3] * x ** 0
12 
13 popt = opt.curve_fit(yfun, x, y, [.01, 1, np.pi/4, .0125])[0]
14 print(popt)

where [.01, 1, np.pi/4, .0125] represents initial guesses for the magnitude a, frequency b, phase c, and average value of the function. The popt variable would have 4 values.

References