Fixing inflexion point estimate using python

2020-06-06 04:55发布

I am trying to find the inflexion points on a curve using python. The data for the curve is here: https://www.dropbox.com/s/rig8frgewde8i5n/fitted.txt?dl=0. Please note that the curve has been fitted to the raw data. Raw data is available here: https://www.dropbox.com/s/1lskykdi1ia1lu7/ww.txt?dl=0

import numpy as np
# Read in array from text file
arr = np.loadtxt(path_to_file)

inflexion_point_1 = np.diff(arr).argmin()
inflexion_point_2 = np.diff(arr).argmax()

These 2 inflexion points are shows as red lines in the attached plot. However, their locations do not seem to be right. The first inflexion point should be close to the area indicated by the black arrow. How do I fix this?

enter image description here

Also, here is a plot of the differential:

plt.axvline(np.gradient(arr[:365]).argmax())

As you can see, the code is behaving as coded i.e. it finds the argmax of np.diff of the array. However, I want to find a position closer to day 110 or so, i.e. about half-way to argmax.

enter image description here

--EDIT--

Also, here is another plot showing thee raw data and the fitted curve (using a quadratic function).

enter image description here

1条回答
狗以群分
2楼-- · 2020-06-06 05:47

Any reason not to use uni-variate spline directly on the gradient?

from scipy.interpolate import UnivariateSpline

#raw data
data = np.genfromtxt('ww.txt')

plt.plot(np.gradient(data), '+')

spl = UnivariateSpline(np.arange(len(data)), np.gradient(data), k=5)
spl.set_smoothing_factor(1000)
plt.plot(spl(np.arange(len(data))), label='Smooth Fct 1e3')
spl.set_smoothing_factor(10000)
plt.plot(spl(np.arange(len(data))), label='Smooth Fct 1e4')
plt.legend(loc='lower left')

max_idx = np.argmax(spl(np.arange(len(data))))
plt.vlines(max_idx, -5, 9, linewidth=5, alpha=0.3)

enter image description here

Also we can solve for the maximum numerically:

In [122]:

import scipy.optimize as so
F = lambda x: -spl(x)
so.fmin(F, 102)
Optimization terminated successfully.
         Current function value: -3.339112
         Iterations: 20
         Function evaluations: 40
Out[122]:
array([ 124.91303558])
查看更多
登录 后发表回答