Piecewise fitting with two linear functions and a breaking point in Python

Question

I want to fit my data using two linear functions (broken power law) with one breaking point which is user given. Currently Im using the curve_fit function from the scipy.optimize module. Here are my datasets frequencies, binned data, errors

Here is my code:

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


freqs=np.loadtxt('binf11.dat')
binys=np.loadtxt('binp11.dat')
errs=np.loadtxt('bine11.dat')


def brkPowLaw(xArray, breakp, slopeA, offsetA, slopeB):
    returnArray = []
    for x in xArray:
        if x <= breakp:
            returnArray.append(slopeA * x + offsetA)
        elif x>breakp:           
            returnArray.append(slopeB * x + offsetA)
    return returnArray

#define initial guesses, breakpoint=-3.2
a_fit,cov=curve_fit(brkPowLaw,freqs,binys,sigma=errs,p0=(-3.2,-2.0,-2.0,-2.0))


modelPredictions = brkPowLaw(freqs, *a_fit) 


plt.errorbar(freqs, binys, yerr=errs, fmt='kp',fillstyle='none',elinewidth=1)
plt.xlim(-5,-2)
plt.plot(freqs,modelPredictions,'r')

The offset of the second linear function is set to be equal to the offset of the first one.

It looks like this works but I get this fit:

Now I thought that the condition in by brkPowLaw function should suffice but it does not. What I want is that the first linear equation is used to fit the data up to a chosen breaking point and then from this breaking point a second linear fit will be done, but without the hump as it shows in the plot because now there it looks like there are two breaking points instead of one and three linear functions for fitting which is not what I expected nor wanted.

What I want is that when the first linear fit ends the second one starts from the point where the first linear fit ended.

I have tried using the numpy.piecewise function with no plausible result, looked into some topics like this or this but I did not manage to make my script work

Thank you for your time

Answer 1

This would be my approach, not with linear but quadratic functions.

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

def soft_step(x, s): ### not my usual np.tanh() ...OMG
    return 1+ 0.5 * s * x / np.sqrt( 1 + ( s * x )**2 )

### for the looks of the data I decided to go for two parabolas with 
### one discontinuity
def fit_func( x, a0, b0, c0, a1, b1, c1, x0, s ):
    out  = ( a0 * x**2 + b0 * x + c0 ) * ( 1 - soft_step( x - x0, s ) )
    out += ( a1 * x**2 + b1 * x + c1 ) * soft_step( x - x0, s )
    return out

### with global parameter for iterative fit
### if using least_squares one could avoid globals
def fit_short( x, a0, b0, c0, a1, b1, c1, x0 ):
    global stepwidth
    return fit_func( x, a0, b0, c0, a1, b1, c1, x0, stepwidth )

### getting data
xl = np.loadtxt( "binf11.dat" )
yl = np.loadtxt( "binp11.dat" )
el = np.loadtxt( "bine11.dat" )

### check for initial values
p0 = [ 0, -2,-11, 0, -2, -9, -3, 10 ]
xth = np.linspace( -5.5, -1.5, 250 )
yth = np.fromiter( ( fit_func(x, *p0 ) for x in xth ), np.float )

### initial fit
sol, pcov = curve_fit( fit_func, xl, yl, sigma=el, p0=p0, absolute_sigma=True )
yft = np.fromiter( ( fit_func( x, *sol ) for x in xth ), np.float )
sol=sol[: -1]

###iterating with fixed and decreasing softness in the step
for stepwidth in range(10,55,5):
    sol, pcov = curve_fit( fit_short, xl, yl, sigma=el, p0=sol, absolute_sigma=True )
    ### printing the step position
    print sol[-1]
yiter = np.fromiter( ( fit_short(x, *sol ) for x in xth ), np.float )
print sol

###plotting
fig = plt.figure()
ax = fig.add_subplot( 1, 1, 1 )
# ~ax.plot( xth, yth ) ### no need to show start parameters
ax.plot( xth, yft ) ### first fit with variable softness
ax.plot( xth, yiter ) ### last fit with fixed softness of 50
ax.errorbar( xl, yl, el, marker='o', ls='' ) ### data
plt.show()

This gives:

-3.1762721614559712
-3.1804393481217477
-3.1822672190583603
-3.183493292415725
-3.1846976088390333
-3.185974760198917
-3.1872472903175266
-3.188427041827035
-3.1894705102541843
[ -0.78797351  -5.33255174 -12.48258537   0.53024954   1.14252783 -4.44589397  -3.18947051]

and

putting the jump at -3.189