如何在python中適合step函數

Question

我有一個關於使用scipy例程（如curve_fit）擬合step函數的問題。我有麻煩使它向量化，例如：

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

xobs=np.linspace(0,10,100) 
yl=np.random.rand(50); yr=np.random.rand(50)+100 
yobs=np.concatenate((yl,yr),axis=0) 

def model(x,rf,T1,T2): 
#1: x=np.vectorize(x) 
    if x<rf: 
     ret= T1 
    else: 
     ret= T2 
    return ret 
#2: model=np.vectorize(model) 
popt, pcov = curve_fit(model, xobs, yobs, [40.,0.,100.]) 

它說

ValueError: The truth value of an array with more than one element is ambiguous. Use a.any() or a.all() 

如果我添加1號或2運行，但並不真正適合的數據：

OptimizeWarning: Covariance of the parameters could not be estimated  category=OptimizeWarning) 
[ 40.   50.51182064 50.51182064] [[ inf inf inf] 
[ inf inf inf] 
[ inf inf inf]] 

有人知道如何解決這個問題嗎？ THX

Answer 1

這就是我所做的。我保留xobs和yobs：

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

xobs=np.linspace(0,10,100) 
yl=np.random.rand(50); yr=np.random.rand(50)+100 
yobs=np.concatenate((yl,yr),axis=0) 

現在，必須生成Heaviside函數。爲了給你這個功能的概述，認爲Heaviside函數的半最大約定：

在Python，這相當於：def f(x): return 0.5 * (np.sign(x) + 1)

樣本情節是：

xval = sorted(np.concatenate([np.linspace(-5,5,100),[0]])) # includes x = 0 
yval = f(xval) 
plt.plot(xval,yval,'ko-') 
plt.ylim(-0.1,1.1) 
plt.xlabel('x',size=18) 
plt.ylabel('H(x)',size=20) 

現在，繪製xobs和yobs給出：

plt.plot(xobs,yobs,'ko-') 
plt.ylim(-10,110) 
plt.xlabel('xobs',size=18) 
plt.ylabel('yobs',size=20) 

注意，比較這兩個圖中，第二曲線是由5個單位和最大偏移增加爲1.0〜100。我推斷爲第二曲線的函數可以是表示如下：

或在Python：(0.5 * (np.sign(x-5) + 1) * 100 = 50 * (np.sign(x-5) + 1)

個結合的田塊產量（其中Fit表示上述擬合函數）

的情節證實我的猜測是正確的。現在，假設您不知道這個正確的擬合函數是如何產生的，則會創建一個廣義擬合函數：def f(x,a,b,c): return a * (np.sign(x-b) + c)，理論上，其中a = 50,b = 5和c = 1。

繼續進行估算：

popt,pcov=curve_fit(f,xobs,yobs,bounds=([49,4.75,0],[50,5,2]))。

現在，bounds = ([lower bound of each parameter (a,b,c)],[upper bound of each parameter])。從技術上講，這意味着49 < a < 50，4.75 < b < 5和0 < c < 2。

下面是我給popt和pcov結果：

pcov表示POPT的估計的協方差。對角線提供了參數估計值[Source]的方差。

結果顯示參數估計值pcov接近理論值。

基本上，廣義Heaviside函數可以表示爲：a * (np.sign(x-b) + c)

這裏是將產生參數估計和對應的協方差的代碼：

import numpy as np 
from scipy.optimize import curve_fit 

xobs = np.linspace(0,10,100) 
yl = np.random.rand(50); yr=np.random.rand(50)+100 
yobs = np.concatenate((yl,yr),axis=0) 

def f(x,a,b,c): return a * (np.sign(x-b) + c) # Heaviside fitting function 

popt, pcov = curve_fit(f,xobs,yobs,bounds=([49,4.75,0],[50,5,2])) 
print 'popt = %s' % popt 
print 'pcov = \n %s' % pcov 

最後，注意的popt的估計和pcov各不相同。

pythonscipy