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

Question

Я использую 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 * covariance(parameters) = sum of errors(residuals) где sigma(parameters) это то, что я хочу.

Как мне связать то, что делает Curve_fit с тем, что я вижу, например? Википедия: http://en.wikipedia.org/wiki/Propagation_of_uncertainty

Answer 1

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

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

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

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

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

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

Answer 2

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

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

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

Answer 3

Математика

Сначала мы начнем с линейной регрессии. Во многих статистических задачах мы предполагаем, что переменные имеют некоторые основные распределения с некоторыми неизвестными параметрами, и оцениваем эти параметры. В линейной регрессии мы предполагаем, что зависимые переменные 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 β).

Я просто выписываю решение, которое

β^ = (Х^Т Х)^-1 Х^Т 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 σ².

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

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

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

σ ^² = Q / (п - р)

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

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

SciPy API

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

В 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 Y.

И

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)^-1 Q / (п - р).

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