Search code examples
pythonnumpycurve-fittingleast-squares

SciPy + Numpy: Finding the slope of a sigmoid curve

I have some data that follow a sigmoid distribution as you can see in the following image: Sigmoid data for 2003

After normalizing and scaling my data, I have adjusted the curve at the bottom using scipy.optimize.curve_fit and some initial parameters:

popt, pcov = curve_fit(sigmoid_function, xdata, ydata, p0 = [0.05, 0.05, 0.05])
>>> print popt
[  2.82019932e+02  -1.90996563e-01   5.00000000e-02]

So popt, according to the documentation, returns *"Optimal values for the parameters so that the sum of the squared error of f(xdata, popt) - ydata is minimized". I understand here that there is no calculation of the slope with curve_fit, because I do not think the slope of this gentle curve is 282, neither is negative.

Then I tried with scipy.optimize.leastsq, because the documentation says it returns "The solution (or the result of the last iteration for an unsuccessful call).", so I thought the slope would be returned. Like this:

p, cov, infodict, mesg, ier = leastsq(residuals, p_guess, args = (nxdata, nydata), full_output=True)
>>> print p
Param(x0=281.73193626250207, y0=-0.012731420027056234, c=1.0069006606656596, k=0.18836680131910222)

But again, I did not get what I expected. curve_fit and leastsq returned almost the same values, with is not surprising I guess, as curve_fit is using an implementation of the least squares method within to find the curve. But no slope back...unless I overlooked something.

So, how to calculate the slope in a point, say, where X = 285 and Y = 0.5?

I am trying to avoid manual methods, like calculating the derivative in, say, (285.5, 0.55) and (284.5, 0.45) and subtract and divide results and so. I would like to know if there is a more automatic method for this.

Thank you all!

EDIT #1

This is my "sigmoid_function", used by curve_fit and leastsq methods:

def sigmoid_function(xdata, x0, k, p0): # p0 not used anymore, only its components (x0, k)
    # This function is called by two different methods: curve_fit and leastsq,
    # this last one through function "residuals". I don't know if it makes sense
    # to use a single function for two (somewhat similar) methods, but there 
    # it goes.

    # p0:
    #   + Is the initial parameter for scipy.optimize.curve_fit. 
    #   + For residuals calculation is left empty
    #   + It is initialized to [0.05, 0.05, 0.05]
    # x0:
    #   + Is the convergence parameter in X-axis and also the shift
    #   + It starts with 0.05 and ends up being around ~282 (days in a year)
    # k:
    #   + Set up either by curve_fit or leastsq
    #   + In least squares it is initially fixed at 0.5 and in curve_fit
    #   + to 0.05. Why? Just did this approach in two different ways and 
    #   + it seems it is working. 
    #   + But honestly, I have no clue on what it represents
    # xdata: 
    #   + Positions in X-axis. In this case from 240 to 365

# Finally I changed those parameters as suggested in the answer. 
# Sigmoid curve has 2 degrees of freedom, therefore, the initial 
# guess only needs to be this size. In this case, p0 = [282, 0.5]


    y = np.exp(-k*(xdata-x0)) / (1 + np.exp(-k*(xdata-x0)))
    return y

def residuals(p_guess, xdata, ydata):
    # For the residuals calculation, there is no need of setting up the initial parameters
    # After fixing the initial guess and sigmoid_function header, remove [] 
    # return ydata - sigmoid_function(xdata, p_guess[0], p_guess[1], [])
    return ydata - sigmoid_function(xdata, p_guess[0], p_guess[1], [])

I am sorry if I made mistakes while describing the parameters or confused technical terms. I am very new with numpy and I have not studied maths for years, so I am catching up again.

So, again, what is your advice to calculate the slope of X = 285, Y = 0.5 (more or less the midpoint) for this dataset? Thanks!!

EDIT #2

Thanks to Oliver W., I updated my code as he suggested and understood a bit better the problem.

There is a final detail I do not fully get. Apparently, curve_fit returns a popt array (x0, k) with the optimum parameters for the fitting:

  • x0 seems to be how shifted is the curve by indicating the central point of the curve
  • k parameter is the slope when y = 0.5, also in the center of the curve (I think!)

Why if the sigmoid function is a growing one, the derivative/slope in popt is negative? Does it make sense?

I used sigmoid_derivative to calculate the slope and, yes, I obtained the same results that popt but with positive sign. 

# Year 2003, 2005, 2007. Slope in midpoint.
k = [-0.1910, -0.2545, -0.2259] # Values coming from popt
slope = [0.1910, 0.2545, 0.2259] # Values coming from sigmoid_derivative function

I know this is being a bit peaky because I could use both. The relevant data is in there but with negative sign, but I was wondering why is this happening.

So, the calculation of the derivative function as you suggested, is only required if I need to know the slope in other points than y = 0.5. Only for midpoint, I can use popt.

Thanks for your help, it saved me a lot of time. :-)

Solution

  • You're never using the parameter p0 you're passing to your sigmoid function. Hence, curve fitting will not have any good measure to find convergence, because it can take any value for this parameter. You should first rewrite your sigmoid function like this:

    def sigmoid_function(xdata, x0, k):

    y = np.exp(-k*(xdata-x0)) / (1 + np.exp(-k*(xdata-x0)))
    return y

    This means your model (the sigmoid) has only two degrees of freedom. This will be returned in popt:

    initial_guess = [282, 1]  # (x0, k): at x0, the sigmoid reaches 50%, k is slope related
popt, pcov = curve_fit(sigmoid_function, xdata, ydata, p0=initial_guess)

    Now popt will be a tuple (or array of 2 values), being the best possible x0 and k.

    To get the slope of this function at any point, to be honest, I would just calculate the derivative symbolically as the sigmoid is not such a hard function. You will end up with:

    def sigmoid_derivative(x, x0, k):
    f = np.exp(-k*(x-x0))
    return -k / f

    If you have the results from your curve fitting stored in popt, you could pass this easily to this function:

    print(sigmoid_derivative(285, *popt))

    which will return for you the derivative at x=285. But, because you ask specifically for the midpoint, so when x==x0 and y==.5, you'll see (from the sigmoid_derivative) that the derivative there is just -k, which can be observed immediately from the curve_fit output you've already obtained. In the output you've shown, that's about 0.19.