Search code examples
python-3.xnumpyscipyderivative

Python - Find min, max, and inflection points of normal curve


I am trying to find the locations (i.e., the x-value) of minimum, start of season, peak growing season, maximum growth, senescence, end of season, minimum (i.e., inflection points) in a vegetation curve. I am using a normal curve here as an example. I did come across few codes to find the change in slope and 1st/2nd order derivative, but not able to implement them for my case. Please direct me if there is any relevant example and your help is appreciated. Thanks!

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

Solution

  • Here is a very simple and not robust method to find the inflection point of a non-noisy curve:

    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();
    

    The actual position of the inflection point is x1 = mean - std for a Gaussian curve.

    For this to work with real data, they have to be smoothed before looking for the max, by applying for example a simple moving average, a gaussian filter or a Savitzky-Golay filter which can directly output the second derivative... the choice of the right filter depends on the data