curve_fit not optimizing one of the parameters

Question

I need to fit with scipy.optimize.curve_fit some data that look like the points in the figure. I use a function y(x) (see def below) which gives a constant y(x)=c for x<x0, otherwise a polynomial (eg a second tilted line y1 = mx+q).

I give a reasonable initial guess for the parameters (x0, c, m, q), as show in the figure. The result from the fit shows that all the parameters are optimized except for the first one x0.

Why so? Is it how I define the function testfit(x, *p), where x0 (=p[0]) appears within another function?

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit

# generate some data:
x = np.linspace(0,100,1000)
y1 = np.repeat(0, 500)
y2 = x[500:] - 50
y = np.concatenate((y1,y2))
y = y + np.random.randn(len(y))


def testfit(x, *p):
    ''' piecewise function used to fit 
        it's a constant (=p[1]) for x < p[0]
        or a polynomial for x > p[0]     
    '''
    x = x.astype(float)
    y = np.piecewise(x, [x < p[0], x >= p[0]], [p[1], lambda x: np.poly1d(p[2:])(x)])
    return y

# initial guess, one horizontal and one tilted line:
p0_guess = (30, 5, 0.3, -10)

popt, pcov = curve_fit(testfit, x, y, p0=p0_guess)

print('params guessed  : '+str(p0_guess))
print('params from fit : '+str(popt))

plt.plot(x,y, '.')
plt.plot(x, testfit(x, *p0_guess), label='initial guess')
plt.plot(x, testfit(x, *popt), label='final fit')
plt.legend()

Output

params guessed  : (30, 5, 0.3, -10) 
params from fit : [ 30. 0.04970411   0.80106256 -34.17194401] 

OptimizeWarning: Covariance of the parameters could not be estimated category=OptimizeWarning)

Answer 1

As suggested by kazemakase, I solved the problem with a smooth transition between the two functions I use to fit (one horizontal line followed by a polynomial). The trick was to multiply one function by sigmoid(x) and the other by 1-sigmoid(x), (where sigmoid(x) is defined below).

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit

x = np.linspace(0,100,1000)
y1 = np.repeat(0, 500)
y2 = x[500:] - 50
y = np.concatenate((y1,y2))
y = y + np.random.randn(len(y))

def testfit(x, *p):
    ''' function to fit the indentation curve 
    p = [x0,c, poly1d_coeffs ]'''
    x = x.astype(float)
    y = p[1]*(1-sigmoid(x-p[0],k=1)) + np.poly1d(p[2:])(x) * sigmoid(x-p[0],k=1)
    return y

def sigmoid(x, k=1):
    return 1/(1+np.exp(-k*x))

p0_guess = (30, 5, 0.3, -10 )
popt, pcov = curve_fit(testfit, x, y, p0=p0_guess)
print('params guessed  : '+str(p0_guess))
print('params from fit : '+str(popt))


plt.figure(1)
plt.clf()
plt.plot(x,y, 'y.')
plt.plot(x, testfit(x, *p0_guess), label='initial guess')
plt.plot(x, testfit(x, *popt), 'k', label='final fit')
plt.legend()

Answer 2

I had a similar problem. I ended up using np.gradient and a convolution to smooth the curve, then plotting it. Something like:

def mov_avg(n, data):
    return np.convolve(data, np.ones((n,))/n, mode='valid')

If you want a more direct approach, you can try this:

def find_change(data):
    def test_flag(pos):
        grad = np.gradient(data) - np.gradient(data).mean()
        return (grad[:pos]<0).sum() + (grad[pos:]>0).sum()
    return np.vectorize(test_flag)(np.arange(len(data)-1)).argmax()

def find_gradient(pos, data):
    return np.gradient(data[:pos]).mean(), np.gradient(data[pos:]).mean()

pos=find_change(x2)
print(pos, find_gradient(pos, data))

The first function calculates the point at which the gradient change by comparing the point gradient against the mean gradient, and finds the point from which the gradients are "mostly positive".

Hope it helps