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

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

Прежде всего, он говорит, что это якобиан, но в примечаниях он также говорит, что "cov_x является якобианским приближением к гессианскому" , так что это на самом деле не якобиан, а гессиан, используя некоторое приближение от якобиана. Какое из этих утверждений верно?

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

Эта матрица должна быть умножена на остаточную дисперсию, чтобы получить ковариацию оценок параметров - см. 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 * ковариация (параметры) = сумма ошибок (остатков), где сигма (параметры) - это то, что я хочу.

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

Ответ 1

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

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

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

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

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

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

Ответ 2

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

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

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