Python：分段多项式曲线以指数拟合

Question

我想复制一个研究结果，其中分布是按以下方式分段拟合的。

1. ln（y_data）= A（x）：拟合中心部分的三阶多项式
2. ln（y_data）= B（x）：拟合的两个尾部的一阶多项式

到目前为止，我使用curve_fit并定义了一个分段函数，并通过要求在边界x0和x1处施加了连续性和平滑性：

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）

我的代码是：

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

当前x0，x1的值是静态的。理想情况下，这些也是自由变量。

满足边界条件：

# 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

屈服

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

如何获得正确的身材？感谢您的帮助！

（经过下采样的）数据由：

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]])

Answer 1

您可以将曲线拟合问题包装到外部优化问题中，在其中优化x0和x1。

尝试：

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

这将产生以下拟合：click here