Search code examples
pythonnumpycurve-fittingpiecewisescipy-optimize

Python: piecewise polynomial curve fit in the exponential

I want to reproduce a research result where the distribution is piecewise fitted the following way.

  1. ln(y_data) = A(x) : 3rd order polynomial for central part of the fit
  2. ln(y_data) = B(x) : 1st order polynomial for both tails of the fit

So far, I used curve_fit and defined a piecewise function and imposed continuity and smoothness at the boundaries x0 and x1 by requiring:

  1. f_1(x0) = f_2(x0) & f_2(x1) = f_3(x1)
  2. f'_1(x0) = f'_2(x0) & f'_2(x1) = f'_3(x1)

My code is: 

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

plt.figure(figsize=(9,6))

def piecewise_func(x, x0, x1, a0, a1, b0, b1, b2, b3, c0, c1):    
    x0 = 0.015
    x1 = 0.06

    b3 = ((c1-(a0+a1*x0-c0-c1*x1)/(x0-x1) + (x0+x1)/(2*(x0-x1)) * (a1-c1) -x1/(x0-x1)*(a1-c1))/(3.*x1**2))/(1.-(-1.5*(x0+x1)**2 + (x0**3-x1**3)/(x0-x1) + 3.*x1*(x0+x1))/(3.*x1**2))
    b2 = (a1 - c1 - 3.*b3*(x0**2 - x1**2))/(2.*(x0-x1))
    b1 = (a0+a1*x0-c0-c1*x1-b2*(x0**2-x1**2)-b3*(x0**3-x1**3))/(x0-x1)
    b0 = a0 + a1*x0 - b1*x0 - b2*x0**2 - b3*x0**3 

    condlist = [x < x0, (x >= x0) & (x < x1), x >= x1]
    funclist = [lambda x: a0 + a1*x, 
                lambda x: b0 + b1*x + b2*x**2 + b3*x**3, 
                lambda x: c0 + c1*x]
    return np.piecewise(x, condlist, funclist)

p , e = optimize.curve_fit(piecewise_func, xvals, np.log(yvals))
x_plot = np.linspace(0., 0.15, len(xvals))
plt.plot(xvals, np.log(yvals), "+")
plt.plot(xvals, (piecewise_func(x_plot, *p)), 'C3-', lw=3)
plt.show()

Currently the values for x0, x1 are static. Ideally, these would be free variables as well.

The boundary conditions are met:

# continuity
print 'continuity f_1(x0)=f_2(x0): ',a0 + a1*x0, '//', b0 + b1*x0 + b2*x0**2. + b3*x0**3. 
print 'continuity f_2(x1)=f_3(x1): ',b0 + b1*x1 + b2*x1**2. + b3*x1**3.,  '//',c0 + c1*x1 ,'\n'
# smoothness for first derivative
print 'smoothness df_1(x0)=df_2(x0): ',a1, '//', b1 + 2.*b2*x0 + 3.*b3*x0**2.
print 'smoothness df_2(x1)=df_3(x1): ',b1 + 2.*b2*x1 + 3.*b3*x1**2.,  '//',c1

yielding

continuity f_1(x0)=f_2(x0):  4.68285214128 // 4.68285214128
continuity f_2(x1)=f_3(x1):  6.53425912688 // 6.53425912688 

smoothness df_1(x0)=df_2(x0):  233.211740443 // 233.211740443
smoothness df_2(x1)=df_3(x1):  -5.51084715909 // -5.51084715909

curve fit MWE How can I obtain a correct fit? Thanks for your help!

The (downsampled) data is given by: 

array([[  2.72580000e-03,   2.00000000e+00],
       [  5.35160000e-03,   6.00000000e+00],
       [  8.35440000e-03,   3.10000000e+01],
       [  1.15380000e-02,   8.10000000e+01],
       [  1.48130000e-02,   1.89000000e+02],
       [  1.81220000e-02,   2.88000000e+02],
       [  2.14320000e-02,   4.31000000e+02],
       [  2.47350000e-02,   5.27000000e+02],
       [  2.80210000e-02,   5.90000000e+02],
       [  3.12730000e-02,   6.59000000e+02],
       [  3.44930000e-02,   7.56000000e+02],
       [  3.76910000e-02,   7.90000000e+02],
       [  4.08610000e-02,   8.36000000e+02],
       [  4.39970000e-02,   8.61000000e+02],
       [  4.71140000e-02,   8.45000000e+02],
       [  5.01930000e-02,   8.38000000e+02],
       [  5.32830000e-02,   8.39000000e+02],
       [  5.63310000e-02,   8.14000000e+02],
       [  5.93720000e-02,   7.61000000e+02],
       [  6.24040000e-02,   7.29000000e+02],
       [  6.54260000e-02,   7.06000000e+02],
       [  6.84390000e-02,   7.05000000e+02],
       [  7.14430000e-02,   6.77000000e+02],
       [  7.44370000e-02,   6.45000000e+02],
       [  7.74210000e-02,   6.82000000e+02],
       [  8.04130000e-02,   6.75000000e+02],
       [  8.34120000e-02,   6.77000000e+02],
       [  8.63850000e-02,   6.63000000e+02],
       [  8.93820000e-02,   6.62000000e+02],
       [  9.23710000e-02,   6.34000000e+02],
       [  9.53670000e-02,   5.87000000e+02],
       [  9.83720000e-02,   5.46000000e+02],
       [  1.01390000e-01,   5.22000000e+02],
       [  1.04390000e-01,   5.10000000e+02],
       [  1.07400000e-01,   4.97000000e+02],
       [  1.10440000e-01,   4.77000000e+02],
       [  1.13480000e-01,   4.82000000e+02],
       [  1.16540000e-01,   4.55000000e+02]])
Solution

  • You could wrap the curve fit problem into an outer optimization problem where you optimize x0 and x1.

    Try:

    def piecewise_func(p,x):
    a0,a1,c0,c1,x0,x1 = p

    b3 = ((c1-(a0+a1*x0-c0-c1*x1)/(x0-x1) + (x0+x1)/(2*(x0-x1)) * (a1-c1) -x1/(x0-x1)*(a1-c1))/(3.*x1**2))/ (1.-(-1.5*(x0+x1)**2 + (x0**3-x1**3)/(x0-x1) + 3.*x1*(x0+x1))/(3.*x1**2))
    b2 = (a1 - c1 - 3.*b3*(x0**2 - x1**2))/(2.*(x0-x1))
    b1 = (a0+a1*x0-c0-c1*x1-b2*(x0**2-x1**2)-b3*(x0**3-x1**3))/(x0-x1)
    b0 = a0 + a1*x0 - b1*x0 - b2*x0**2 - b3*x0**3

    condlist = [x < x0, (x >= x0) & (x < x1), x >= x1]
    funclist = [lambda x: a0 + a1*x,
                lambda x: b0 + b1*x + b2*x**2 + b3*x**3,
                lambda x: c0 + c1*x]
    return np.piecewise(x, condlist, funclist)


def optFunc(p,extra):
    x,y = extra
    return piecewise_func(p,x)-y


fit,ier = optimize.leastsq(optFunc,[0.,1.,0.,1.,0.01,0.1],[x,np.log(y)])
print fit

plt.plot(x,np.log(y),'k+')
plt.plot(x,piecewise_func(fit,x))

    This should yield the following fit: click here