Prawidłowe dopasowanie za pomocą scipy curve_fit z błędami x?

Question

Próbuję dopasować histogram z niektórych danych w nim przy użyciu scipy.optimize.curve_fit. Jeśli chcę dodać błąd w y, mogę to zrobić, nakładając na dopasowanie weight. Ale jak zastosować błąd w x (np. Błąd związany z binowaniem w przypadku histogramów)?

Moje pytanie dotyczy również błędów w x podczas regresji liniowej z curve_fit lub polyfit; Wiem, jak dodawać błędy w y, ale nie w x.

Oto przykład (częściowo z 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--') 

Teraz, to pasuje (gdy nie powiedzie) uważa Y błędy sigma_n, ale nie znalazłem sposób zmusić go do rozważenia sigma_x. Przeskanowałem kilka wątków na liście dyskusyjnej scipy i dowiedziałem się, jak używać wartości absolute_sigma i posta na Stackoverflow o asymmetrical errors, ale nic o błędach w obu kierunkach. Czy można to osiągnąć?

Answer 1

scipy.optmize.curve_fit wykorzystuje standardową nieliniową optymalizację najmniejszych kwadratów, a zatem minimalizuje tylko odchylenie w zmiennych odpowiedzi. Jeśli chcesz uwzględnić błąd w zmiennej niezależnej, możesz wypróbować scipy.odr, która wykorzystuje regresję odległości ortogonalnej. Jak sama nazwa wskazuje minimalizuje zarówno w zmiennych niezależnych, jak i zależnych.

Zobacz przykład poniżej. Parametr fit_type określa, czy scipy.odr ma pełną ODR (fit_type=0) lub optymalizację najmniejszych kwadratów (fit_type=2).

EDIT

Chociaż przykład pracował to nie ma większego sensu, ponieważ dane y została obliczona na hałaśliwych x danych, które właśnie spowodowało w nierównomiernie rozmieszczone zmiennej indepenent. Zaktualizowałem próbkę, która teraz pokazuje również, jak używać RealData, która pozwala określić standardowy błąd danych zamiast wag.

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)