Как рассчитать доверительный интервал для наименьших квадратов (scipy.optimize.leastsq) в Python?
доверительный интервал с минимальным соответствием в scipy python
Ответы (3)
Я бы использовал метод начальной загрузки.
См. Здесь: http://phe.rockefeller.edu/LogletLab/whitepaper/node17.html
Простой пример для шумного гауссовского языка:
x = arange(-10, 10, 0.01)
# model function
def f(p):
mu, s = p
return exp(-(x-mu)**2/(2*s**2))
# create error function for dataset
def fff(d):
def ff(p):
return d-f(p)
return ff
# create noisy dataset from model
def noisy_data(p):
return f(p)+normal(0,0.1,len(x))
# fit dataset to model with least squares
def fit(d):
ff = fff(d)
p = leastsq(ff,[0,1])[0]
return p
# bootstrap estimation
def bootstrap(d):
p0 = fit(d)
residuals = f(p0)-d
s_residuals = std(residuals)
ps = []
for i in range(1000):
new_d = d+normal(0,s_residuals,len(d))
ps.append(fit(new_d))
ps = array(ps)
mean_params = mean(ps,0)
std_params = std(ps,0)
return mean_params, std_params
data = noisy_data([0.5, 2.1])
mean_params, std_params = bootstrap(data)
print "95% confidence interval:"
print "mu: ", mean_params[0], " +/- ", std_params[0]*1.95996
print "sigma: ", mean_params[1], " +/- ", std_params[1]*1.95996
Я не уверен, что вы имеете в виду под доверительным интервалом.
В общем, leastsq
мало что знает о функции, которую вы пытаетесь минимизировать, поэтому он не может дать доверительный интервал. Однако он возвращает оценку гессиана, другими словами, обобщение 2-х производных на многомерные проблемы.
Как указано в строке документации функции, вы можете использовать эту информацию вместе с остатками (разницей между вашим подобранным решением и фактическими данными) для вычисления ковариации оценок параметров, которая является локальным предположением доверительного интервала.
Обратите внимание, что это только локальная информация, и я подозреваю, что, строго говоря, вы можете прийти к выводу, только если ваша целевая функция строго выпуклая. Никаких подтверждений или ссылок на это утверждение у меня нет :).
Самый простой способ оценки доверительного интервала (ДИ) - умножить стандартные ошибки (стандартное отклонение) на константу. Чтобы вычислить константу, вам необходимо знать количество степеней свободы (DOF) и уровень достоверности, для которого вы хотите рассчитать CI. Полученный таким образом КИ иногда называют асимптотическим КИ. Подробнее об этом можно прочитать в разделе «Подбор моделей к биологическим данным с использованием линейной и нелинейной регрессии» Мотульски и Кристопулоса (книги Google). Эта же книга (или очень похожая) доступна бесплатно в качестве руководства для авторского программного обеспечения.
Вы также можете прочитать как рассчитать CI с помощью библиотеки C ++ Boost.Math. В этом примере CI рассчитывается для распределения одной переменной. В случае аппроксимации методом наименьших квадратов глубина резкости не N -1, а N-M, где M - количество параметров. В Python должно быть легко сделать то же самое.
Это простейшая оценка. Я не знаю метода начальной загрузки, предложенного zephyr, но он может быть более надежным, чем метод, о котором я писал.