Python 中的二阶导数 - scipy/numpy/pandas

Posted 2023-02-23

【中文标题】Python 中的二阶导数 - scipy/numpy/pandas

【英文标题】：Second Derivative in Python - scipy/numpy/pandas

【发布时间】：2017-03-06 16:47:09

【问题描述】：

我正在尝试在 python 中使用两个 numpy 数据数组进行二阶导数。

例如，有问题的数组如下所示：

import numpy as np

x = np.array([ 120. ,  121.5,  122. ,  122.5,  123. ,  123.5,  124. ,  124.5,
        125. ,  125.5,  126. ,  126.5,  127. ,  127.5,  128. ,  128.5,
        129. ,  129.5,  130. ,  130.5,  131. ,  131.5,  132. ,  132.5,
        133. ,  133.5,  134. ,  134.5,  135. ,  135.5,  136. ,  136.5,
        137. ,  137.5,  138. ,  138.5,  139. ,  139.5,  140. ,  140.5,
        141. ,  141.5,  142. ,  142.5,  143. ,  143.5,  144. ,  144.5,
        145. ,  145.5,  146. ,  146.5,  147. ])

y = np.array([  1.25750000e+01,   1.10750000e+01,   1.05750000e+01,
         1.00750000e+01,   9.57500000e+00,   9.07500000e+00,
         8.57500000e+00,   8.07500000e+00,   7.57500000e+00,
         7.07500000e+00,   6.57500000e+00,   6.07500000e+00,
         5.57500000e+00,   5.07500000e+00,   4.57500000e+00,
         4.07500000e+00,   3.57500000e+00,   3.07500000e+00,
         2.60500000e+00,   2.14500000e+00,   1.71000000e+00,
         1.30500000e+00,   9.55000000e-01,   6.65000000e-01,
         4.35000000e-01,   2.70000000e-01,   1.55000000e-01,
         9.00000000e-02,   5.00000000e-02,   2.50000000e-02,
         1.50000000e-02,   1.00000000e-02,   1.00000000e-02,
         1.00000000e-02,   1.00000000e-02,   1.00000000e-02,
         1.00000000e-02,   1.00000000e-02,   5.00000000e-03,
         5.00000000e-03,   5.00000000e-03,   5.00000000e-03,
         5.00000000e-03,   5.00000000e-03,   5.00000000e-03,
         5.00000000e-03,   5.00000000e-03,   5.00000000e-03,
         5.00000000e-03,   5.00000000e-03,   5.00000000e-03,
         5.00000000e-03,   5.00000000e-03])

我目前有f(x) = y，我想要d^2 y / dx^2。

在数值上，我知道我可以对函数进行插值并分析求导，也可以使用higher order finite-differences。我认为有足够的数据可以使用，如果一个或另一个被认为更快、更准确等。

我查看了 np.interp() 和 scipy.interpolate 没有成功，因为这返回了一个拟合（线性或三次）样条曲线，但不知道如何获得该点的导数。

非常感谢任何指导。

【问题讨论】：

你看过np.diff吗？

我担心我的数据点间距不均匀。

【参考方案1】：

您可以使用 scipy 的一维 Splines 函数对数据进行插值。计算样条有一个方便的derivative 方法来计算导数。

对于您示例的数据，使用 UnivariateSpline 给出以下拟合

import matplotlib.pyplot as plt
from scipy.interpolate import UnivariateSpline

y_spl = UnivariateSpline(x,y,s=0,k=4)

plt.semilogy(x,y,'ro',label = 'data')
x_range = np.linspace(x[0],x[-1],1000)
plt.semilogy(x_range,y_spl(x_range))

合身看起来相当不错，至少在视觉上是这样。您可能想试验UnivariateSpline 使用的参数。

样条拟合的二阶导数可以简单地获得为

y_spl_2d = y_spl.derivative(n=2)

plt.plot(x_range,y_spl_2d(x_range))

结果看起来有些不自然（如果您的数据对应于某个物理过程）。您可能想要更改样条拟合参数、改进数据（例如，提供更多样本、执行较少噪声的测量），或者决定使用分析函数对数据进行建模并执行曲线拟合（例如，使用 sicpy 的 curve_fit )

【讨论】：

这个数据应该代表一个概率密度函数。规范化这条曲线和应用一些规则（比如没有负值）等的最佳方法是什么？

我认为没有标准答案，因为使用通用插值方法的方法在施加约束方面的选择有限。原则上，您需要从头开始制定和解决约束优化问题。您可能希望首先标准化您的数据，因为y_spl.integral(x[0],x[-1]) 大约是 80，这当然不是 pdf 的有效值。

这个答案与两次使用np.diff 的比率有何不同？在数值上是好是坏？

【参考方案2】：

通过有限差分，数组上 x 的每个平均值的 y 的一阶导数由下式给出：

dy=np.diff(y,1)
dx=np.diff(x,1)
yfirst=dy/dx

而x对应的值为：

xfirst=0.5*(x[:-1]+x[1:])

对于第二个订单，再次执行相同的过程：

dyfirst=np.diff(yfirst,1)
dxfirst=np.diff(xfirst,1)
ysecond=dyfirst/dxfirst

xsecond=0.5*(xfirst[:-1]+xfirst[1:])

【讨论】：

np.diff(np.diff([x*x for x in range(0,10)])) = [2,2,2..]