Python - 寻找正态分布曲线的最小值、最大值和拐点

Question

我正在尝试找到植被曲线上最小值、季节开始、生长高峰、最大生长、衰老、季节结束和最小值（即拐点）的位置（即x值）。这里以正常曲线为例。我找到了一些代码来查找斜率变化和一阶/二阶导数，但无法将它们用于我的情况。如果有任何相关示例，请指引我，非常感谢您的帮助。谢谢！

## Version 2 code
import matplotlib.pyplot as plt
import numpy as np
from  scipy.stats import norm

x_min = 0.0
x_max = 16.0

mean = 8
std = 2

x = np.linspace(x_min, x_max, 100)
y = norm.pdf(x, mean, std)

# Slice the group in 3
def group_in_threes(slicable):
    for i in range(len(slicable)-2):
        yield slicable[i:i+3]

# Locate the change in slope
def turns(L):
    for index, three in enumerate(group_in_threes(L)):
        if (three[0] > three[1] < three[2]) or (three[0] < three[1] > three[2]):
            yield index + 1

# 1st inflection point estimation
dy = np.diff(y, n=1) # first derivative
idx_max_dy = np.argmax(dy)
ix = list(turns(dy))
print(ix)

# All inflection point estimation
dy2 = np.diff(dy, n=2) # Second derivative?
idx_max_dy2 = np.argmax(dy2)
ix2 = list(turns(dy2))
print(ix2)

# Graph
plt.plot(x, y)
#plt.plot(x[ix], y[ix], 'or', label='estimated inflection point')
plt.plot(x[ix2], y[ix2], 'or', label='estimated inflection point - 2')
plt.xlabel('x'); plt.ylabel('y'); plt.legend(loc='best');

- pyPN

---

我认为解决方案将取决于实际数据，例如噪声的级别类型...如果曲线看起来像正态分布曲线，为什么不拟合一个正态分布曲线（使用scipy.optimize.curve_fit），然后从获得的均值和方差计算所需的特征？这些特征，即季节开始、生长高峰、最大生长、衰老、季节结束等，如何被“数学”定义？ - xdze2

感谢@xdze2！我使用的曲线只是代表性的。我主要是在寻找在曲线中先于拐点的顶点列表。一种方法是使用二阶导数，并查找从“+ve”到“-ve”或反之亦然的符号变化。我只是不知道如何做。 - pyPN

Answer 1

这是一种非常简单但不太健壮的方法，可用于查找非嘈杂曲线的拐点：

import matplotlib.pyplot as plt
import numpy as np
from  scipy.stats import norm

x_min = 0.0
x_max = 16.0

mean = 8
std = 2

x = np.linspace(x_min, x_max, 100)
y = norm.pdf(x, mean, std)

# 1st inflection point estimation
dy = np.diff(y) # first derivative
idx_max_dy = np.argmax(dy)

# Graph
plt.plot(x, y)
plt.plot(x[idx_max_dy], y[idx_max_dy], 'or', label='estimated inflection point')
plt.xlabel('x'); plt.ylabel('y'); plt.legend();

实际拐点位置为高斯曲线的 x1 = mean - std 。

为了处理真实数据，必须在查找最大值之前对其进行平滑处理，例如应用简单移动平均、高斯滤波器或Savitzky-Golay 滤波器，它们可以直接输出二次导数...选择正确的滤波器取决于数据。

- xdze2

---

感谢@xdze2！很高兴知道Savitzky-Golay滤波器输出二阶导数，并且正如你所说，滤波取决于数据。这是我正在寻找的东西（上面的版本2代码），不确定我是否正确处理了二阶导数部分？ - pyPN