В Scipy как и почему curve_fit вычисляет ковариацию оценок параметров

Question

Я использую scipy.optimize.leastsq для подгонки некоторых данных. Я хотел бы получить некоторые доверительные интервалы по этим оценкам, поэтому я смотрю на вывод cov_x, но документация очень неясна в отношении того, что это такое и как получить из этого ковариационную матрицу для моих параметров.

Прежде всего, это говорит, что это якобиан, но в notes он также говорит, что «cov_x является якобиан приближение к Гессу», так что это на самом деле не якобиан но Hessian используя некоторую аппроксимацию от якобиану , Какое из этих утверждений верно?

Во-вторых, это предложение для меня это сбивает с толку:

Эта матрица должна быть умножена на остаточной дисперсии, чтобы получить ковариации оценок параметров - см curve_fit.

Я действительно посмотрите на исходный код для curve_fit, где они делают:

s_sq = (func(popt, *args)**2).sum()/(len(ydata)-len(p0)) 
pcov = pcov * s_sq 

, который соответствует умножению cov_x на s_sq, но я не могу найти это уравнение в какой-либо ссылки. Может ли кто-нибудь объяснить, почему это уравнение верно? Моя интуиция говорит мне, что это должно быть наоборот, поскольку cov_x должен быть производным (якобиан или гессиан), поэтому я думал: cov_x * covariance(parameters) = sum of errors(residuals), где sigma(parameters) - это то, что я хочу.

Как подключить вещь curve_fit делает то, что я вижу, например. wikipedia: http://en.wikipedia.org/wiki/Propagation_of_uncertainty#Non-linear_combinations

Answer 1

ОК, я думаю, что нашел ответ. Сначала решение: cov_x * s_sq - это просто ковариация параметров, которые вы хотите. Взятие sqrt диагональных элементов даст вам стандартное отклонение (но будьте осторожны с ковариациями!).

Остаточная дисперсия = приведенная цифра = s_sq = sum [(f (x) -y)^2]/(N-n), где N - количество точек данных, а n - количество параметров подгонки. Reduced chi square.

Причина моей путаницы в том, что cov_x, заданный наименьшим значением, на самом деле не является тем, что называется cov (x) в других местах, а является уменьшенным cov (x) или дробным cov (x). Причина, по которой он не обнаруживается ни в одной из других ссылок, заключается в том, что это простое масштабирование, которое полезно в числовых вычислениях, но не относится к учебнику.

О Гессиан против Якобиана, документация плохо сформулирована. В обоих случаях вычисляется гессиан, что очевидно, так как якобиан равен нулю в минимуме. Они имеют в виду, что они используют приближение к якобиану, чтобы найти гессиан.

Еще одно примечание. Кажется, что результат curve_fit фактически не учитывает абсолютный размер ошибок, а учитывает только относительный размер предоставленных сигм. Это означает, что возвращенные pcov не изменяются, даже если ошибки будут меняться на миллион. Это, конечно, не так, но, похоже, это стандартная практика. Matlab делает то же самое при использовании набора инструментов Curve.Правильная процедура описана здесь: https://en.wikipedia.org/wiki/Linear_least_squares_(mathematics)#Parameter_errors_and_correlation

Это кажется довольно простым для этого, как только оптимальный найден, по крайней мере, для линейных наименьших квадратов.

Answer 2

Я нашел это решение во время поиска аналогичного вопроса, и у меня есть лишь небольшое улучшение в ответе Хансхархоффа. Полный вывод из lesssq предоставляет обратное значение infodict, которое содержит infodict ['fvec'] = f (x) -y. Таким образом, чтобы вычислить уменьшенные х-квадрат = (в приведенном выше обозначений)

s_sq = (infodict['fvec']**2).sum()/ (N-n)

КСТАТИ. Спасибо HansHarhoff за большую часть тяжелого подъема, чтобы решить эту проблему.