如何在Python中應用分段線性擬合？

Question

我試圖適應分段線性擬合，如圖1所示爲一個數據組

由該位線設定獲得這個數字。我試圖使用代碼應用分段線性擬合：

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


x = np.array([1, 2, 3, 4, 5, 6, 7, 8, 9, 10 ,11, 12, 13, 14, 15]) 
y = np.array([5, 7, 9, 11, 13, 15, 28.92, 42.81, 56.7, 70.59, 84.47, 98.36, 112.25, 126.14, 140.03]) 


def linear_fit(x, a, b): 
    return a * x + b 
fit_a, fit_b = optimize.curve_fit(linear_fit, x[0:5], y[0:5])[0] 
y_fit = fit_a * x[0:7] + fit_b 
fit_a, fit_b = optimize.curve_fit(linear_fit, x[6:14], y[6:14])[0] 
y_fit = np.append(y_fit, fit_a * x[6:14] + fit_b) 


figure = plt.figure(figsize=(5.15, 5.15)) 
figure.clf() 
plot = plt.subplot(111) 
ax1 = plt.gca() 
plot.plot(x, y, linestyle = '', linewidth = 0.25, markeredgecolor='none', marker = 'o', label = r'\textit{y_a}') 
plot.plot(x, y_fit, linestyle = ':', linewidth = 0.25, markeredgecolor='none', marker = '', label = r'\textit{y_b}') 
plot.set_ylabel('Y', labelpad = 6) 
plot.set_xlabel('X', labelpad = 6) 
figure.savefig('test.pdf', box_inches='tight') 
plt.close()  

但是，這給了我適合圖。 2，我試圖玩弄價值觀，但沒有改變，我無法得到適合的上限。對我來說最重要的要求是如何讓Python獲得漸變點。本質上，我希望Python能夠識別並在適當的範圍內適合兩種線性擬合。這怎麼可以在Python中完成？

Answer 1

您可以使用numpy.piecewise()創建分段函數，然後使用curve_fit()，這裏是代碼

from scipy import optimize 
import matplotlib.pyplot as plt 
import numpy as np 
%matplotlib inline 

x = np.array([1, 2, 3, 4, 5, 6, 7, 8, 9, 10 ,11, 12, 13, 14, 15], dtype=float) 
y = np.array([5, 7, 9, 11, 13, 15, 28.92, 42.81, 56.7, 70.59, 84.47, 98.36, 112.25, 126.14, 140.03]) 

def piecewise_linear(x, x0, y0, k1, k2): 
    return np.piecewise(x, [x < x0], [lambda x:k1*x + y0-k1*x0, lambda x:k2*x + y0-k2*x0]) 

p , e = optimize.curve_fit(piecewise_linear, x, y) 
xd = np.linspace(0, 15, 100) 
plt.plot(x, y, "o") 
plt.plot(xd, piecewise_linear(xd, *p)) 

輸出：

Answer 2

你可以做一個spline interpolation方案來執行分段線性插值並找出曲線的轉折點。的二階導數將在轉折點最高（對於單調增加的曲線），並且可以與順序> 2.

import numpy as np 
import matplotlib.pyplot as plt 
from scipy import interpolate 

x = np.array([1, 2, 3, 4, 5, 6, 7, 8, 9, 10 ,11, 12, 13, 14, 15]) 
y = np.array([5, 7, 9, 11, 13, 15, 28.92, 42.81, 56.7, 70.59, 84.47, 98.36, 112.25, 126.14, 140.03]) 

tck = interpolate.splrep(x, y, k=2, s=0) 
xnew = np.linspace(0, 15) 

fig, axes = plt.subplots(3) 

axes[0].plot(x, y, 'x', label = 'data') 
axes[0].plot(xnew, interpolate.splev(xnew, tck, der=0), label = 'Fit') 
axes[1].plot(x, interpolate.splev(x, tck, der=1), label = '1st dev') 
dev_2 = interpolate.splev(x, tck, der=2) 
axes[2].plot(x, dev_2, label = '2st dev') 

turning_point_mask = dev_2 == np.amax(dev_2) 
axes[2].plot(x[turning_point_mask], dev_2[turning_point_mask],'rx', 
      label = 'Turning point') 
for ax in axes: 
    ax.legend(loc = 'best') 

plt.show() 

Answer 3

使用numpy.interp一個樣條內插，它返回一個來計算維分段線性插值到具有離散數據點給定值的函數。

Answer 4

擴展@ binoy-pilakkat的回答。

你應該使用numpy.interp：

import numpy as np 
import matplotlib.pyplot as plt 

x = np.array(range(1,16), dtype=float) 
y = np.array([5, 7, 9, 11, 13, 15, 28.92, 
      42.81, 56.7, 70.59, 84.47, 
      98.36, 112.25, 126.14, 140.03], dtype=float) 

yinterp = np.interp(x, x, y) # simple as that 

plt.plot(x, y, 'bo') 
plt.plot(x, yinterp, 'g-') 
plt.show() 

Answer 5

具有線性表達做外推法和內插一切。

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

x = np.array([1, 2, 3, 4, 5, 6, 7, 8, 9, 10 ,11, 12, 13, 14, 15], dtype=float) 
y = np.array([5, 7, 9, 11, 19, 15, 28.92, 42.81, 56.7, 79, 84.47, 98.36, 112.25, 126.14, 140.03]) 

def piecewise_linear(x, x0, y0): 
    y = np.zeros(len(x)) 
    def get_slope(xi, xj, yi, yj): 
     return (yj - yi)/(xj - xi) 
    def get_intercept(m, xi, yi): 
     return yi - m*xi 
    def left_piecewise_linear(x, x0, y0): 
     m_left = get_slope(x0[0], x0[1], y0[0], y0[1]) 
     b_left = get_intercept(m_left, x0[0], y0[0]) 
     return m_left * x + b_left 
    def right_piecewise_linear(x, x0, y0): 
     m_right = get_slope(x0[-1], x0[-2], y0[-1], y0[-2]) 
     b_right = get_intercept(m_right, x0[-1], y0[-1]) 
     return m_right * x + b_right 
    for i in range(0,len(x)): 
     if x[i] < x0[0]: 
      y[i] = left_piecewise_linear(x[i], x0, y0) 
     elif x[i] > x0[-1]: 
      y[i] = right_piecewise_linear(x[i], x0, y0) 
     else: 
      y[i] = np.interp(x[i], x0, y0) 
    return y   

xd = np.linspace(-15, 25, 100)  
plt.plot(x, y, "o") 
plt.plot(xd, piecewise_linear(xd, x, y)) 
plt.savefig('img.png')