Правильная установка с scipy curve_fit, включая ошибки в x?

Я пытаюсь установить гистограмму с некоторыми данными в ней с помощью 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 о асимметричных ошибках, но ничего об ошибках в в обоих направлениях. Можно ли достичь?

Ответы

Ответ 1

scipy.optmize.curve_fit использует стандартную нелинейную оптимизацию наименьших квадратов и, следовательно, только минимизирует отклонение в переменных ответа. Если вы хотите, чтобы ошибка в независимой переменной считалась, вы можете попробовать scipy.odr, которая использует ортогональную дистанционную регрессию. Как видно из его названия, он минимизирует как независимые, так и зависимые переменные.

Взгляните на образец ниже. Параметр fit_type определяет, выполняет ли scipy.odr полную ODR (fit_type=0) или оптимизацию наименьших квадратов (fit_type=2).

ИЗМЕНИТЬ

Хотя пример работал, это не имело большого смысла, так как данные y были рассчитаны на шумные данные x, что приводило к неравномерно разнесенной независимой переменной. Я обновил образец, который теперь также показывает, как использовать 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)

fit to noisy data