在Python中擬合分段函數

Question

我試圖將一個分段定義的函數擬合到Python中的數據集中。我已經搜索了很長一段時間，但我還沒有找到答案，可能與否。

要了解我正在嘗試做什麼，請看下面的示例（這不適用於我）。在這裏，我試圖將一個移位的絕對值函數（f（x）= | x-p |）擬合到以p爲擬合參數的數據集中。

import scipy.optimize as so 
import numpy as np 

def fitfunc(x,p): 
    if x>p: 
     return x-p 
    else: 
     return -(x-p) 

fitfunc = np.vectorize(fitfunc) #vectorize so you can use func with array 

x=np.arange(1,10) 
y=fitfunc(x,6)+0.1*np.random.randn(len(x)) 

popt, pcov = so.curve_fit(fitfunc, x, y) #fitting routine that gives error 

有什麼辦法可以在Python中完成這項工作嗎？

R中這樣做的一個方法是：

# Fit of a absolute value function f(x)=|x-p| 

f.lr <- function(x,p) { 
    ifelse(x>p, x-p,-(x-p)) 
} 
x <- seq(0,10) # 
y <- f.lr(x,6) + rnorm (length(x),0,2) 
plot(y ~ x) 
fit.lr <- nls(y ~ f.lr(x,p), start = list(p = 0), trace = T, control = list(warnOnly = T,minFactor = 1/2048)) 
summary(fit.lr) 
coefficients(fit.lr) 
p.fit <- coefficients(fit.lr)["p"] 
x_fine <- seq(0,10,length.out=1000) 
lines(x_fine,f.lr(x_fine,p.fit),type='l',col='red') 
lines(x,f.lr(x,6),type='l',col='blue') 

後甚至更多的研究，我發現這樣做的一種方式。在這個解決方案中，我不喜歡我必須自己定義錯誤函數的事實。此外，我不太確定它爲什麼必須採用這種lambda風格。因此，任何建議或更復雜的解決方案都非常受歡迎。

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

def fitfunc(p,x): return x - p if x > p else p - x 

def array_fitfunc(p,x): 
    y = np.zeros(x.shape) 
    for i in range(len(y)): 
     y[i]=fitfunc(x[i],p) 
    return y 

errfunc = lambda p, x, y: array_fitfunc(p, x) - y # Distance to the target function 

x=np.arange(1,10) 
x_fine=np.arange(1,10,0.1) 
y=array_fitfunc(6,x)+1*np.random.randn(len(x)) #data with noise 

p1, success = so.leastsq(errfunc, -100, args=(x, y), epsfcn=1.) # -100 is the initial value for p; epsfcn sets the step width 

plt.plot(x,y,'o') # fit data 
plt.plot(x_fine,array_fitfunc(6,x_fine),'r-') #original function 
plt.plot(x_fine,array_fitfunc(p1[0],x_fine),'b-') #fitted version 
plt.show() 

Answer 1

要在這裏完成這件事，我會分享我自己的最終解決問題。爲了保持接近我原來的問題，你只需要自己定義向量化函數，而不是使用np.vectorize。

import scipy.optimize as so 
import numpy as np 

def fitfunc(x,p): 
    if x>p: 
     return x-p 
    else: 
     return -(x-p) 

fitfunc_vec = np.vectorize(fitfunc) #vectorize so you can use func with array 

def fitfunc_vec_self(x,p): 
    y = np.zeros(x.shape) 
    for i in range(len(y)): 
    y[i]=fitfunc(x[i],p) 
    return y 


x=np.arange(1,10) 
y=fitfunc_vec_self(x,6)+0.1*np.random.randn(len(x)) 

popt, pcov = so.curve_fit(fitfunc_vec_self, x, y) #fitting routine that gives error 
print popt 
print pcov 

輸出：

[ 6.03608994] 
[[ 0.00124934]] 

Answer 2

你不能簡單地用

def fitfunc2(x, p): 
    return np.abs(x-p) 

然後產生使用類似

>>> x = np.arange(1,10) 
>>> y = fitfunc2(x,6) + 0.1*np.random.randn(len(x)) 
>>> 
>>> so.curve_fit(fitfunc2, x, y) 
(array([ 5.98273313]), array([[ 0.00101859]])) 

開關功能和/或積木一樣where更換樹枝代替fitfunc，這應該擴展到更復雜的表達式，而無需致電vectorize。

[PS：您最小二乘示例中的errfunc不需要是lambda。你可以寫

def errfunc(p, x, y): 
    return array_fitfunc(p, x) - y 

相反，如果你喜歡。]