scipy curve_fit的正確擬合，包括x中的錯誤？

Question

我試圖用scipy.optimize.curve_fit來擬合一些數據的直方圖。如果我想在y中添加一個錯誤，我可以簡單地通過將weight應用於適合。但是如何應用x中的錯誤（即在直方圖情況下由於分箱而導致的錯誤）？

我的問題也適用於x的錯誤，當用curve_fit或polyfit進行線性迴歸時;我知道如何在y中添加錯誤，但不在x中。

下面的例子（部分來自matplotlib documentation）：

import numpy as np 
import pylab as P 
from scipy.optimize import curve_fit 

# create the data histogram 
mu, sigma = 200, 25 
x = mu + sigma*P.randn(10000) 

# define fit function 
def gauss(x, *p): 
    A, mu, sigma = p 
    return A*np.exp(-(x-mu)**2/(2*sigma**2)) 

# the histogram of the data 
n, bins, patches = P.hist(x, 50, histtype='step') 
sigma_n = np.sqrt(n) # Adding Poisson errors in y 
bin_centres = (bins[:-1] + bins[1:])/2 
sigma_x = (bins[1] - bins[0])/np.sqrt(12) # Binning error in x 
P.setp(patches, 'facecolor', 'g', 'alpha', 0.75) 

# fitting and plotting 
p0 = [700, 200, 25] 
popt, pcov = curve_fit(gauss, bin_centres, n, p0=p0, sigma=sigma_n, absolute_sigma=True) 
x = np.arange(100, 300, 0.5) 
fit = gauss(x, *popt) 
P.plot(x, fit, 'r--') 

現在，這種配合（當它不失敗）不會考慮在y錯誤sigma_n，但我還沒有找到一種方法，使其考慮sigma_x。我掃描了scipy郵件列表上的一些線程，並發現瞭如何使用absolute_sigma值和Stackoverflow上關於asymmetrical errors的帖子，但沒有涉及雙向錯誤。是否有可能實現？

Answer 1

scipy.optmize.curve_fit使用標準的非線性最小二乘法優化，因此只能最小化響應變量的偏差。如果你想要考慮自變量的錯誤，你可以試試scipy.odr，它使用正交距離迴歸。顧名思義，它將獨立和因變量都降到最低。

看看下面的示例。 fit_type參數確定scipy.odr是否執行完整的ODR（fit_type=0）或最小二乘法優化（fit_type=2）。

編輯

儘管示例的工作並沒有多大意義，因爲Y數據是在嘈雜的X數據，剛剛產生了一個不等距indepenent變量計算。我更新了樣本，該樣本現在還顯示如何使用RealData，它允許指定數據的標準誤而不是權重。

from scipy.odr import ODR, Model, Data, RealData 
import numpy as np 
from pylab import * 

def func(beta, x): 
    y = beta[0]+beta[1]*x+beta[2]*x**3 
    return y 

#generate data 
x = np.linspace(-3,2,100) 
y = func([-2.3,7.0,-4.0], x) 

# add some noise 
x += np.random.normal(scale=0.3, size=100) 
y += np.random.normal(scale=0.1, size=100) 

data = RealData(x, y, 0.3, 0.1) 
model = Model(func) 

odr = ODR(data, model, [1,0,0]) 
odr.set_job(fit_type=2) 
output = odr.run() 

xn = np.linspace(-3,2,50) 
yn = func(output.beta, xn) 
hold(True) 
plot(x,y,'ro') 
plot(xn,yn,'k-',label='leastsq') 
odr.set_job(fit_type=0) 
output = odr.run() 
yn = func(output.beta, xn) 
plot(xn,yn,'g-',label='odr') 
legend(loc=0)