Python: вычисление доверительных интервалов для нелинейного подбора

Question

У меня есть некоторые данные x и y data -> Просто в качестве примера

x_data = [5, 6, 6.2, 6.3, 6.7]
y_data = [5.2, 6.8, 8.2, , 8.9] 

, и теперь я хочу добавить нелинейную функцию

def func(x, a):
     return 0.2* a**2 * x**2

со scipy. Для этого я использую

prams, params_cov = optimize.curve_fit(func, x_data, y_data)

Но у меня также есть полосы ошибок для данных x и y. Так что я не хочу помещать только строку, а полосу ошибок с fillbetween.

x_err = [0.2, 0.1, 0.4, 0.5]
y_err = [0.1, 0.2, 0.6, 0.03]

Но для этого мне нужны доверительные интервалы относительно полос ошибок. Кто-нибудь знает, как их вычислить? Или как подобрать функцию, относящуюся к полосам ошибок и получить из нее ошибку для подгонки?

С наилучшими пожеланиями

Answer 1

Использование оценки соответствия параметров по методу Монте-Карло для учета ваших ошибок будет выглядеть следующим образом:

import numpy as np
from scipy.optimize import curve_fit

def myfunc(x, a):
     return 0.2* a**2 * x**2

x_data = np.array([5, 6.2, 6.3, 6.7])
y_data = np.array([5.2, 6.8, 8.2, 8.9])
x_err = np.array([0.2, 0.1, 0.4, 0.5])
y_err = np.array([0.1, 0.2, 0.6, 0.03])

trials = 1000
my_guess = (2)

fit_params = np.zeros((1,trials))

for i in range(trials):
    x_sim = x_data + np.random.normal(scale=x_err)
    y_sim = y_data + np.random.normal(scale=y_err)

    params, params_cov = curve_fit(myfunc, x_sim, y_sim, my_guess)

    fit_params[0,i] = params

mu = np.mean(fit_params)
ci = np.std(fit_params) * 1.96 / np.sqrt(np.size(x_data))
print ("Mean a = {:.3f} +/- {:.3f} (95% CI)".format(mu, ci))

Выход: Mean a = 0.987 +/- 0.031 (95% CI)

Примечание. Предполагается, что ваши ошибки распределяются в соответствии с нормальным распределением.

Чаще всего используют t-критерий с учетом степеней свободы и дисперсии для вашего подобранного параметра, который можно рассчитать следующим образом:

dof = np.size(x_data) - 1 # degrees of freedom:
# calculate student-t value
a = 0.05 #(1-0.95, 95% CI)
tval = t.ppf(1.0-a/2, dof)

params, params_cov = curve_fit(myfunc, x_data, y_data, 1)

ci = tval*np.sqrt(params_cov)[0]
print ("Mean a = {:.3f} +/- {:.3f} (95% CI)".format(params.item(), ci.item()))

Выход: Mean a = 0.990 +/- 0.056 (95% CI)