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.

 

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out /  Change )

Google photo

You are commenting using your Google account. Log Out /  Change )

Twitter picture

You are commenting using your Twitter account. Log Out /  Change )

Facebook photo

You are commenting using your Facebook account. Log Out /  Change )

Connecting to %s