## 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 * angle ** 3) + (beta * angle ** 2) + (beta * angle) + beta torque_fit.append(y) elif order == 4: for angle in angle_fit: # calculate torques from 4th order polynomial y = (beta * angle ** 4) + (beta * angle ** 3) + (beta * angle ** 2) + (beta * angle) + beta 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.

• philipperemy89

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

```np.apply_along_axis(func, 0, angle_fit, a, b, c)
```

Like

• Joanna Diong

Thanks Philippe for the helpful improvements. Appreciate the interest and feedback.

Like

• philipperemy89

Nice article!

One thing to improve on:

``````if mathexp == 'poly':
# perform least squares polynomial fit
beta = np.polyfit(angle, torque, order)
torque_fit = np.polyval(beta, angle_fit)
``````

Will be much faster and elegant than naively looping.

Thanks!

Like

• Nice article