0
votes

I have a set of data, y is angular orientation, and x is the timestamp for each point of y.

The entire data set has many segments for angular orientation. Inorder to do curve fitting, I have split the data into their respective segments, storing each segment as an numpy array.

I then find a polynomial fit using numpy.polyfit to find a curve fit for each segment of the data. However because my data is purely experimental, I have no knowledge of which polynomial degree to use for numpy.polyfit, thus I iterate through a range of polynomial degrees till I find the highest polynomial degree possible.

This is my code:

# Experimental data stored in lists: time_aralist and orienrad_aralist
# List stores the segments as arrays

fig = plt.figure()

# Finding curve fit
fittime_aralist, fitorienrad_aralist, fitorienrad_funclist = [], [], []
for j in range(len(time_aralist)):
    z, res, _, _, _ = np.polyfit(time_aralist[j], orienrad_aralist[j], 200, full=True)
    orienrad_func = np.poly1d(z)
    fittime_ara = np.linspace(time_aralist[j][0], time_aralist[j][-1], 10000)
    fitorienrad_ara = orienrad_func(fittime_ara)

    # Appending to list
    fittime_aralist.append(fittime_ara)
    fitorienrad_aralist.append(fitorienrad_ara)
    fitorienrad_funclist.append(orienrad_func)

# Plotting experimental data
for j in range(len(time_aralist)):
    plt.plot(time_aralist[j], orienrad_aralist[j], 'ro')
for j in range(len(fittime_aralist)):
    plt.plot(fittime_aralist[j], fitorienrad_aralist, 'k-')

This is what my plot looks like. Centered in the plot is one segment.

The black lines give the attempted curve fit, and the red dots are the experimental points. enter image description here

As seen in the diagram, the black curve fitted lines do not really fit the data points very well. My region of interest is the middle region of the segment, however the region is not fitted well as well, despite using the highest polynomial degree possible.

Could anyone provide me with any supplementary techniques, or an alternative code that could fit the data better?

1
The question is, do you have a model to fit to the data? Or do you just want to produce a nice line through the points? I.e. do you need a fit? If so, one would need to know which possible fitting functions to consider. If not, one shouldn't use fitting at all.ImportanceOfBeingErnest
I require a curve fit as I am interested in the first derivative of the data. I have no knowledge of the model to fit.Tian
My point is that if you don't have a model, there is nothing to fit. Is having a dense first derivative the only aim here?ImportanceOfBeingErnest
Yes that is the main aim.Tian

1 Answers

1
votes

Cubic interpolation

What about a cubic interpolation of the data?

import numpy as np
from scipy.interpolate import interp1d
import matplotlib.pyplot as plt

x = np.linspace(6, 13, num=40)
y = 3 + 2.*x+np.sin(x)/2.+np.sin(4*x)/3.+ np.random.rand(len(x))/3.

fig, ax = plt.subplots()
ax.scatter(x,y, s=5, c="crimson")

f = interp1d(x, y, kind='cubic')
xdens = np.linspace(6, 13, num=400)
ydens = f(xdens)

ax.plot(xdens, ydens, label="interpolation")
ax.legend()
ax2 = ax.twinx()
yderiv =  np.diff(ydens)/np.diff(xdens)
ax2.plot(xdens[:-1],yderiv, color="C2", label="derivative")

ax2.legend()
plt.show()

enter image description here

Spline interpolation

The same can be achieved with a Spline interpolation. The advantage would be that scipy.interpolate.splrep allows to use the parameter s to smoothen the result and it also allows to evaluate directly the derivative of the spline.

import numpy as np
import scipy.interpolate
import matplotlib.pyplot as plt

x = np.linspace(6, 13, num=40)
y = 3 + 2.*x+np.sin(x)/2.+np.sin(4*x)/3.+ np.random.rand(len(x))/3.

fig, ax = plt.subplots()
ax.scatter(x,y, s=5, c="crimson")

tck = scipy.interpolate.splrep(x, y, s=0.4)
xdens = np.linspace(6, 13, num=400)
ydens = scipy.interpolate.splev(xdens, tck, der=0)

ax.plot(xdens, ydens, label="spline interpolation")
ax.legend(loc=2)

ax2 = ax.twinx()
yderiv = scipy.interpolate.splev(xdens, tck, der=1)
ax2.plot(xdens, yderiv, color="C2", label="derivative")

ax2.legend(loc=4)
plt.show()

enter image description here