Сравнение curve_fit и scipy.odr - абсолютная сигма - PullRequest
Новая
       106

Сравнение curve_fit и scipy.odr - абсолютная сигма

1 голос
/ 19 июня 2020

В настоящее время я хочу сопоставить данные с ошибками в x и y и im, используя пакет scipy.odr, чтобы получить свои результаты. Мне просто интересно, как правильно использовать ошибки sx и sy. Вот пример.

Предположим, я измеряю напряжение V для разных токов I. Итак, я измерил 1 В с погрешностью + - 0,1 В и так далее. Поэтому, если я предполагаю, что в текущем измерении нет ошибки, я мог бы использовать spipy.curve_fit следующим образом. Absolute_sigma установлено в True, потому что моя абсолютная ошибка составляет + - 0,1 В.

Итак, я получаю: 

from scipy.optimize import curve_fit
y=[1,2.5,3,4]
x=[1,2,3,4]
yerr=[0.1,0.1,0.1,0.1]

def func(x, a, b):
    return a*x+b

popt, pcov = curve_fit(func, x, y,sigma=yerr,absolute_sigma=True)

print(popt)
print(np.sqrt(np.diag(pcov)))

[0.95 0.25]
[0.04472136 0.12247449]

На втором этапе я хочу использовать odr-package с ошибками в как ток, так и напряжение. Согласно документации это должно использоваться следующим образом: sx и sy - это ошибки для моих данных измерений. Поэтому я должен предположить, что получаю аналогичные результаты с curve_fit, если я использую очень маленькую ошибку для sx. 

from scipy.odr import *
x_err = [0.0000000001]*x.__len__()
y_err = yerr

def linear_func(p, x):
   m, c = p
   return m*x+c

linear_model = Model(linear_func)

data = RealData(x, y, sx=x_err, sy=y_err)

odr = ODR(data, linear_model, beta0=[0.4, 0.4])

out = odr.run()

out.pprint()

Beta: [0.94999996 0.24999994]
Beta Std Error: [0.13228754 0.36228459]

Но, как вы можете видеть, результат отличается от curve_fit выше с absolute_sigma = True.

Использование тех же данных с curve_fit и absolute_sigma = False приводит к тем же результатам, что и ODR-Fit. 

popt, pcov = curve_fit(func, x, y,sigma=yerr,absolute_sigma=False)

print(popt)
print(np.sqrt(np.diag(pcov)))

[0.95 0.25]
[0.13228756 0.36228442]

Итак, я предполагаю, что ODR-Fit на самом деле не заботится о моих абсолютных ошибках, как это делают curve_fit и absolute_sigma = True. Есть ли способ сделать это или я что-то упускаю?

1 Ответ

3 голосов
/ 08 июля 2020

Параметр absolute_sigma=True в curve_fit() выдает матрицу реальной ковариации в том смысле, что np.sqrt(np.diag(pcov)) дает стандартное отклонение в 1 сигму как ошибку, как определено, например, в числовых рецептах, т. Е. , это может быть единица измерения, например метр или около того. Очень полезная сводка содержит kmpfit

ODR, однако выдает стандартную ошибку, полученную из масштабированной ковариационной матрицы. См. Здесь: Как вычислить стандартную ошибку из результатов ODR? или пример ниже.

Это масштабирование с помощью ODR выполняется таким образом, что сокращенный chi2 - вычисляется с входными весами, масштабируемыми в таким же образом - дает приблизительно 1. Или из curve_fit docs :

Эта константа задается требованием, чтобы сокращенный chisq для оптимальных параметров появлялся при использовании масштабированной сигмы, равной единице. .

Остается вопрос: что на самом деле означает sd_beta? Это стандартная ошибка , см., Например, Стандартная ошибка среднего значения в сравнении со стандартным отклонением (могут существовать условия, при которых величина обоих одинакова. См. Комментарии ниже) см. Также предыдущее обсуждение здесь: https://mail.python.org/pipermail/scipy-user/2013-February/034196.html

Теперь, посредством масштабирования pcov с уменьшенным chi2, для ошибок параметров получается тот же результат a) curve_fit(..., absolute_sigma=False) и b) ODR, которые априори относительные ошибки: 

# calculate chi2
residuals = (y - func(x, *popt))
chi_arr =  residuals / yerr
chi2_red = (chi_arr**2).sum() / (len(x)-len(popt))
print('red. chi2:\t\t\t\t', chi2_red)
print('np.sqrt(np.diag(pcov) * chi2_red):\t', np.sqrt(np.diag(pcov) * chi2_red))

, что дает: 

red. chi2:                           8.75
np.sqrt(np.diag(pcov) * chi2_red):   [0.13228756 0.36228442]

Обратите внимание, что, однако, для curve_fit(..., absolute_sigma=True) и ODR ковариационные матрицы pcov из curve_fit и cov_beta из вывода ODR остаются такими же. Здесь становится очевидным недостаток обратного и обратного масштабирования!

относительная ошибка теперь, конечно, подразумевает, что если входные ошибки масштабируются, величина относительных ошибок остается: 

# scaled uncertainty
yerr = np.asfarray([0.1, 0.1, 0.1, 0.1]) * 2

popt, pcov = curve_fit(func, x, y, sigma=yerr, absolute_sigma=False)

print('popt:\t\t\t', popt)
print('np.sqrt(np.diag(pcov)):\t', np.sqrt(np.diag(pcov)))

с то же вывод: 

popt:                    [0.95 0.25]
np.sqrt(np.diag(pcov)):  [0.13228756 0.36228442]
...