Fitting polynomial or exponential curves to biological data in Python
Raw data are not always pretty. Data can have different patterns, be noisy, and vary from trial to trial. We usually collect data to measure, or more accurately, estimate an effect or phenomenon. At times, it is necessary to fit a mathematical expression to raw data in order to estimate the underlying effect or phenomenon. We can then use this mathematical expression to determine outcomes of interest. How is this done?
In this tutorial, we use the file “torque_angle.csv” which is a dataset of passive ankle torque and angle data, available for download here. Passive joint torque usually has a positive exponential relationship to passive angle: angle increases as torque increases. At low torques, a one-unit increase in torque increases joint angle substantially. At high torques, a one-unit increase in torque only produces a small increase in joint angle.
We would like to fit a 4th order polynomial or an exponential function to the raw data and for each type of function, read off the angle when torque is 10 Nm, and plot the raw and fitted data. Here is a function to fit either type of mathematical expression to the data, using the sampling rate freq.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 |
import numpy as np
from scipy.optimize import curve_fit
import pandas as pd
import matplotlib.pyplot as plt
def fit(freq, angle, torque, torque_val, mathexp='poly', order=4):
# find min, max values in angle array
angle_min, angle_max = np.min(angle), np.max(angle)
# calculate fitted angle values based on sampling rate
angle_fit = np.arange(angle_min, angle_max, 1 / freq)
torque_fit = []
if mathexp == 'poly':
# perform least squares polynomial fit
beta = np.polyfit(angle, torque, order)
if order == 3:
for angle in angle_fit:
# calculate torques from 3rd order polynomial
y = (beta[0] * angle ** 3) + (beta[1] * angle ** 2) + (beta[2] * angle) + beta[3]
torque_fit.append(y)
elif order == 4:
for angle in angle_fit:
# calculate torques from 4th order polynomial
y = (beta[0] * angle ** 4) + (beta[1] * angle ** 3) + (beta[2] * angle ** 2) + (beta[3] * angle) + beta[4]
torque_fit.append(y)
torque_fit = np.array(torque_fit)
elif mathexp == 'exp':
# define exponential function
def func(x, a, b, c):
return a * np.exp(-b * x) + c
# initialise starting values: these take some fiddling around
p0 = 5, 0.01, -10
# fit exponential function to data, return optimal values and estimated covariance for parameters
popt, pcov = curve_fit(func, angle, torque, p0, maxfev=10000)
a, b, c = popt
# calculate predicted torque using coefficients and angle
for angle in angle_fit:
y = a * np.exp(-b * angle) + c
torque_fit.append(y)
torque_fit = np.array(torque_fit)
# return indices to sort array, then sort angle and torque
sort_order = angle_fit.argsort()
# find predicted angle at known torque
angle_fit = angle_fit[sort_order]
torque_fit = torque_fit[sort_order]
angle_val = angle_fit[torque_fit.searchsorted(torque_val, 'left')]
# return fitted torque and angle, and predicted angle
return angle_fit, torque_fit, angle_val
|
When the mathematical expression (i.e. mathexp) is specified as polynomial (line 13), we can fit either 3rd or 4th order polynomials to the data, but 4th order is the default (line 7). We use the np.polyfit function to fit a polynomial curve to the data using least squares (line 19 or 24).
Fitting exponential curves is a little trickier. When the mathematical expression is specified as exponential (line 27), the exponential function first has to be defined and its parameters returned. How this function is expressed depends on the physiological relationship of the data; e.g. if the data have a positive exponential form, then the function has to be written as positive exponential. Next, since exponential functions can vary substantially, it is helpful to define and initialise starting values so the fits are better. Finally, we use Scipy’s curve_fit function to fit the exponential curve to the data using non-linear least squares (line 34).
The fitted angle and torque arrays are sorted, then ankle angle at 10 Nm is read off from the fitted curves. Now, we read the data into a Pandas dataframe and perform the fits.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 |
# read data, initialise variables
df = pd.read_csv('torque_angle.csv')
freq = 2000
torque_val = 10
# fit torque-angle data with 4th order polynomial, plot
angle, torque, angle_val = fit(freq, df.angle, df.torque, torque_val, mathexp='poly', order=4)
print('Angle at 10 Nm (polynomial):', angle_val)
fig = plt.figure()
ax = fig.add_subplot(111)
plt.plot(df.angle, df.torque, 'k', label='raw data')
plt.plot(angle, torque, 'b', label='fitted curve')
plt.ylabel('Torque (Nm)')
plt.xlabel('Angle (deg)')
plt.legend()
plt.savefig('fig1.png', dpi=300)
# fit torque-angle data with exponential, plot
angle, torque, angle_val = fit(freq, df.angle, df.torque, torque_val, mathexp='exp')
print('Angle at 10 Nm (exponential):', angle_val)
fig = plt.figure()
ax = fig.add_subplot(111)
plt.plot(df.angle, df.torque, 'k', label='raw data')
plt.plot(angle, torque, 'b', label='fitted curve')
plt.ylabel('Torque (Nm)')
plt.xlabel('Angle (deg)')
plt.legend()
plt.savefig('fig2.png', dpi=300)
|
It is always important to visually check how well the fitted curves match the raw data. Here are the data fitted using a 4th order polynomial: ankle angle at 10 Nm is 93.34 deg.
Figure 1:
Here are the data fitted using an exponential curve: ankle angle at 10 Nm is 93.30 deg.
Figure 2:
Both types of functions fit the data pretty well, and the predicted angles are identical to 1 decimal place. Who would have thought math and Python could be so handy!
Summary
We fitted polynomial and exponential functions to ankle torque-angle data, and used the fitted curves to read off angle at a given torque.
Scientifically Sound 
Same for the exponential part,
you can replace
for angle in angle_fit: y = a * np.exp(-b * angle) + c torque_fit.append(y) torque_fit = np.array(torque_fit)by
LikeLike
Thanks Philippe for the helpful improvements. Appreciate the interest and feedback.
LikeLike
Nice article!
One thing to improve on:
Will be much faster and elegant than naively looping.
Thanks!
LikeLike
Nice article
Instead of doing
y = (beta[0] * angle ** 4) + (beta[1] * angle ** 3) + (beta[2] * angle ** 2) + (beta[3] * angle) + beta[
4]
You can simply do np.polyval(beta, angle_fit)
Values will be the same and I bet it will be much faster.
LikeLike