Python - найти минимальные, максимальные и точки перегиба нормальной кривой

Я пытаюсь найти места (т. е. значение x) минимума, начала сезона, пика вегетационного периода, максимального роста, старения, конца сезона, минимума (т. е. точки перегиба) на кривой вегетации. В качестве примера я использую нормальную кривую. Я столкнулся с несколькими кодами, чтобы найти изменение наклона и производной 1-го и 2-го порядка, но не смог реализовать их для своего случая. Пожалуйста, направьте меня, если есть какой-либо соответствующий пример, и ваша помощь приветствуется. Спасибо!

## 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');

python-3.x numpy scipy derivative
person pyPN    schedule 27.08.2019    source источник
comment
Я думаю, что решение будет зависеть от реальных данных, например, от уровня типа шума... Если кривые выглядят как гауссовы, почему бы не подобрать гауссову (используя scipy.optimize.curve_fit), а затем вычислить нужные функции из полученного среднего значения и дисперсии? Как математически определяются признаки, т. е. начало сезона, пик вегетационного периода, максимальный рост, старение, конец сезона...?   -  person xdze2    schedule 28.08.2019
comment
Спасибо @xdze2! Кривая, которую я использую, является просто репрезентативной. В основном я ищу список вершин, которые предшествуют точкам перегиба кривой. Один из способов — использовать вторую производную и искать изменение знака с +ve на -ve или наоборот. Я просто не знаю, как это сделать.   -  person pyPN    schedule 28.08.2019

Ответы (1)


arrow_upward
1
arrow_downward

Вот очень простой и ненадежный метод найти точку перегиба кривой без шума:

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 для кривой Гаусса.

Чтобы это работало с реальными данными, их необходимо сгладить перед поиском максимума, применив, например, простое скользящее среднее, фильтр Гаусса или фильтр Савицкого-Голея, который может напрямую выводить вторую производную... выбор правильного фильтра зависит от данных

person xdze2    schedule 28.08.2019
comment
Спасибо @xdze2! Приятно знать, что Савицки- Фильтр Голея выводит вторую производную, и, как вы сказали, фильтрация зависит от данных. Вот то, что я ищу (код версии 2 выше), не уверен, правильно ли я сделал вторую производную часть? - person pyPN; 28.08.2019