Lắp đúng với đường cong scipy bao gồm các lỗi trong x?

Question

Tôi đang cố gắng để phù hợp với một biểu đồ với một số dữ liệu trong nó bằng cách sử scipy.optimize.curve_fit. Nếu tôi muốn thêm lỗi trong y, tôi có thể chỉ cần làm như vậy bằng cách áp dụng weight cho phù hợp. Nhưng làm thế nào để áp dụng các lỗi trong x (i. E lỗi do binning trong trường hợp của biểu đồ)?

Câu hỏi của tôi cũng áp dụng cho các lỗi trong x khi thực hiện hồi quy tuyến tính với curve_fit hoặc polyfit; Tôi biết cách thêm lỗi trong y, nhưng không phải trong x.

Dưới đây là một ví dụ (một phần từ 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--') 

Bây giờ, phù hợp này (khi nó không thất bại) không xem xét các y-lỗi sigma_n, nhưng tôi đã không tìm thấy một cách để làm cho nó xem xét sigma_x. Tôi quét một vài chủ đề trên mailing list scipy và phát hiện ra làm thế nào để sử dụng giá trị absolute_sigma và một bài đăng trên Stackoverflow về asymmetrical errors, nhưng không về lỗi trong cả hai hướng. Có thể đạt được không?

Answer 1

scipy.optmize.curve_fit sử dụng tiêu chuẩn phi tuyến tính phương nhỏ nhất tối ưu hóa và do đó chỉ giảm thiểu độ lệch trong các biến phản ứng. Nếu bạn muốn có một lỗi trong biến độc lập được coi là bạn có thể thử scipy.odr sử dụng hồi quy khoảng cách trực giao. Như tên gọi của nó cho thấy nó giảm thiểu trong cả hai biến độc lập và phụ thuộc.

Hãy nhìn vào mẫu dưới đây. Tham số fit_type sẽ xác định xem scipy.odr không ODR đầy đủ (fit_type=0) hoặc ít nhất là hình vuông tối ưu hóa (fit_type=2).

EDIT

Mặc dù ví dụ làm việc nó không có ý nghĩa nhiều, vì các dữ liệu y đã được tính toán trên ồn ào x dữ liệu, mà chỉ dẫn đến một biến indepenent không đều cách nhau. Tôi cập nhật các mẫu mà bây giờ cũng cho thấy làm thế nào để sử dụng RealData cho phép xác định lỗi tiêu chuẩn của dữ liệu thay vì trọng lượng.

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)