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

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

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

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

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

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

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

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

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

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

Математика

Сначала мы начнем с линейной регрессии. Во многих статистических задачах мы предполагаем, что переменные имеют некоторые базовые распределения с некоторыми неизвестными параметрами, и мы оцениваем эти параметры. В линейной регрессии мы предполагаем, что зависимые переменные y_i имеют линейную связь с независимыми переменными x_ij:

y_i = x_i1β₁ + ... + x_ipβ_p > + σε_i, i = 1, ..., n.

где ε_i имеет независимое стандартное нормальное распределение, β_j - это p неизвестных параметров, а также неизвестно σ. Мы можем записать это в матричной форме:

Y = X β + σε,

где Y, β и ε — вектор-столбец. Чтобы найти наилучшее β, мы минимизируем сумму квадратов

S = (Y - X β)^T (Y - X β).

Я просто пишу решение, которое

β^ = (X^T X)^-1 X^T Y.

Если мы видим Y как конкретные наблюдаемые данные, β ^ является оценкой β при этом наблюдении. С другой стороны, если мы рассматриваем Y как случайную величину, оценка β^ также становится случайной величиной. Таким образом, мы можем увидеть, какова ковариация β^.

Поскольку Y имеет многомерное нормальное распределение, а β^ является линейным преобразованием Y, β^ также имеет многомерное нормальное распределение. Ковариационная матрица β^ равна

Cov(β^) = (X^T X)^-1 X^T Cov(Y) ((X^T > X)^-1 X^T)^T = (X ^T X)^{-1 σ2.}

Но здесь σ неизвестен, поэтому его тоже нужно оценить. Если мы позволим

Q = (Y - X β^)^T (Y - X β^),

можно доказать, что Q / σ² имеет распределение хи-квадрат с n - p степенями свободы (при этом Q не зависит от β^). Это делает

σ^² = Q / (n - p)

несмещенная оценка σ². Таким образом, окончательная оценка Cov(β^) равна

(X^T X)^-1Q/(n - p).

SciPy-API

curve_fit наиболее удобен, второе возвращаемое значение pcov — это всего лишь оценка ковариации β^, то есть окончательный результат (X^T X)^-1 Q/ (н - р) выше.

В leastsq второе возвращаемое значение cov_x равно (X^T X)^-1. Из выражения S мы видим, что X^T X является гессианом S (точнее, половиной гессиана), поэтому в документе говорится, что cov_x является обратным гессиану. Чтобы получить ковариацию, нужно умножить cov_x на Q/(n - p).

Нелинейная регрессия

В нелинейной регрессии y_i зависит от параметров нелинейно:

y_i = f(x_i, β₁, ..., β_{p) + σεi.}

Мы можем вычислить частные производные от f по β_j, поэтому она становится приблизительно линейной. Тогда расчет в основном такой же, как линейная регрессия, за исключением того, что нам нужно итеративно аппроксимировать минимум. На практике алгоритм может быть более сложным, например, алгоритм Левенберга-Марквардта, который по умолчанию имеет значение curve_fit.

Подробнее о предоставлении Sigma

Этот раздел посвящен параметрам sigma и absolute_sigma в curve_fit. Для базового использования curve_fit, когда у вас нет предварительных знаний о ковариации Y, вы можете игнорировать этот раздел.

Абсолютная сигма

В приведенной выше линейной регрессии дисперсия y_i равна σ и неизвестна. Если вы знаете разницу. Вы можете предоставить его curve_fit через параметр sigma и установить absolute_sigma=True.

Предположим, что предоставленная вами матрица sigma равна Σ. то есть

Y ~ N(X β, Σ).

Y имеет многомерное нормальное распределение со средним X β и ковариацией Σ. Мы хотим максимизировать вероятность Y. Из функции плотности вероятности Y, которая эквивалентна минимизации

S = (Y - X β)^T Σ^-1 (Y - X β).

Решение

β^ = (X^T Σ^-1 X)^-1 X^T Σ^{-1< /sup> Ю.}

И

Cov(β^) = (X^T Σ^-1 X)^-1.

Приведенные выше β^ и Cov(β^) являются возвращаемыми значениями curve_fit с absolute_sigma=True.

Относительная сигма

В некоторых случаях вы не знаете точную дисперсию y_i, но знаете относительную связь между различными y_i, например дисперсию y ₂ в 4 раза превышает дисперсию y₁. Затем вы можете передать sigma и установить absolute_sigma=False.

Этот раз

Y ~ N(X β, Σσ)

с заданной известной матрицей Σ и неизвестным числом σ. Целевая функция, которую нужно минимизировать, такая же, как абсолютная сигма, поскольку σ является константой, и, следовательно, оценщик β^ тот же. Но ковариация

Cov(β^) = (X^T Σ^-1 X)^-1 σ²,

содержит неизвестное σ. Для оценки σ пусть

Q = (Y - X β^)^T Σ^-1 (Y - X β^).

Опять же, Q / σ² имеет распределение хи-квадрат с n - p степенями свободы.

Оценка Cov(β^) равна

(X^T Σ^-1 X)^-1Q/(n - p).

И это второе возвращаемое значение curve_fit с absolute_sigma=False.