Руководство

Предположим, что для -го наблюдения связь между независимой переменной и зависимой переменной выражается следующей формулой:

где  — нелинейная функция модели, которая зависит от независимой переменной и вектора параметров . Чтобы найти параметр , который лучше всего соответствует нашим данным, мы выбираем параметр так, чтобы свести к минимуму сумму квадратов остатков в данных, то есть решаем следующую задачу:

Учитывая, что функция нелинейная, задача нахождения оптимального значения не имеет аналитического решения. Для решения задачи наименьших квадратов необходимо использовать вычислительные инструменты, а именно LsqFit.jl в данном руководстве.

Одним из примеров нелинейной модели является экспоненциальная модель, которая принимает одноэлементную предикторную переменную . Функция модели имеет следующий вид:

Модель в этом случае описывается так:

Для подгонки данных с помощью LsqFit.jl передайте определенную функцию модели (m), данные (tdata и ydata) и начальное значение параметра (p0) в curve_fit(). В настоящее время LsqFit.jl поддерживает только алгоритм Левенберга-Марквардта.

Он возвращает составной тип LsqFitResult с рядом любопытных значений:

Якобиан векторной функции определяется как матрица со следующими элементами:

Таким образом, матрица имеет вид:

Якобиан экспоненциальной функции модели относительно :

По умолчанию с целью аппроксимации якобиана для алгоритма подгонки данных и вычисления ковариации используется метод конечных разностей (дополнительные сведения см. в описании модуля NLSolversBase.jl). Чтобы получить результаты быстрее и, возможно, с большей точностью, можно также передать функцию, вычисляющую якобиан, в curve_fit().

Нелинейную функцию можно аппроксимировать как линейную путем разложения в ряд Тейлора:

где  — это фиксированный вектор,  — вектор очень небольшой величины, а  — градиент при .

Рассмотрим векторную функцию остатков

с элементами:

Линейную аппроксимацию каждого элемента в таком случае можно записать в следующем виде:

Так как -я строка равна результату транспонирования градиента , векторную функцию можно аппроксимировать так:

Это линейная функция от , так как  — это фиксированный вектор.

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

Рассмотрим решение задачи наименьших квадратов , которое представляет собой локальный минимизатор нелинейной задачи:

Зададим как фиксированную точку в линейной аппроксимации: и . Вектор параметров в окрестности можно выразить как . Локальная аппроксимация для задачи наименьших квадратов имеет следующий вид:

Это, по сути, линейная задача наименьших квадратов:

где , и . Решим уравнение при частных производных, равных . Аналитическое решение выглядит так:

Ковариационная матрица для аналитического решения:

Обратите внимание, что  — это вектор остатков в точке наилучшего соответствия с элементами . Точка очень близка к , поэтому ее можно заменить .

Если допустить, что погрешности в каждой выборке являются независимыми и имеют нормальное распределение с нулевым средним значением и одинаковой дисперсией, то есть , ковариационная матрица линейной аппроксимации имеет следующий вид:

где можно оценить как сумму остатков квадратов, деленную на количество степеней свободы:

В LsqFit.jl для вычисления ковариационной матрицы используется QR-разложение как обладающее большей вычислительной устойчивостью. Оно имеет следующую форму:

vcov() вычисляет ковариационную матрицу соответствия:

Тогда среднеквадратичная погрешность равна корню квадратному из каждого диагонального элемента ковариационной матрицы. stderror() возвращает среднеквадратичную погрешность каждого параметра.

margin_error() вычисляет произведение среднеквадратичной погрешности и критического значения каждого параметра на определенном уровне значимости (по умолчанию 5 %) из t-распределения. Допустимый предел погрешности на уровне значимости 10 % можно вычислить так:

confint() возвращает доверительный интервал каждого параметра на определенном уровне значимости, то есть, по сути, оценочное значение ± допустимый предел погрешности. Чтобы получить доверительный интервал на уровне значимости 10 %, выполните следующий код:

curve_fit() также принимает весовой параметр (wt) для выполнения метода взвешенных наименьших квадратов и общих наименьших квадратов, где параметр минимизирует взвешенную сумму остатков квадратов.

Весовой параметр (wt) представляет собой массив или матрицу весов для каждой выборки. Для выполнения метода взвешенных наименьших квадратов передайте массив весов [w_1, w_2, ..., w_n] или матрицу весов W:

Задача нахождения взвешенных наименьших квадратов принимает вид:

Или в матричной форме:

где

 — векторная функция остатков со следующими элементами:

Алгоритм LsqFit.jl выдает решение задачи наименьших квадратов .

Если вес представляет собой вектор, то в LsqFit.jl функция остатков, передаваемая в levenberg_marquardt(), имеет другой формат:

Если вес задан матрицей, производится разложение Холецкого, то есть, по сути, из матрицы извлекается квадратный корень:

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

Если задать и , линейная аппроксимация задачи нахождения взвешенных наименьших квадратов принимает вид:

Аналитическое решение линейной аппроксимации:

Допустим, что погрешности в каждой выборке являются независимыми и имеют нормальное распределение с нулевым средним значением и разной дисперсией (гетероскедастическую погрешность), то есть , где:

Мы знаем дисперсию погрешности и задаем вес как величину, обратную дисперсии (оптимальный вес), то есть :

Ковариационная матрица теперь будет иметь следующий вид.

Если бы было известно только отношение различных дисперсий, то есть , ковариационная матрица имела бы такую форму:

где является оценочным. Если в этом случае задать , результат будет таким же, как и в версии без веса. Однако curve_fit() в настоящее время не поддерживает такую реализацию. curve_fit() предполагает, что вес является матрицей, обратной ковариационной матрице погрешностей, а не отношением ковариационной матрицы погрешностей, то есть ковариация оценочного параметра вычисляется как covar = inv(J'*fit.wt*J).

Передайте вектор 1 ./ var(ε) или матрицу inv(covar(ε)) в качестве весового параметра (wt) в функцию curve_fit():

Допустим, что погрешности в каждой выборке являются скоррелированными и имеют нормальное распределение с нулевым средним значением и разной дисперсией (гетероскедастическую и автокоррелированную погрешность), то есть .

Если задать матрицу весов как обратную ковариационной матрице погрешностей (оптимальный вес), то есть , то мы получим ковариационную матрицу параметров:

Передайте матрицу inv(covar(ε)) в качестве весового параметра (wt) в функцию curve_fit():

В большинстве случаев дисперсии погрешностей неизвестны. Для выполнения метода взвешенных наименьших квадратов необходимо сначала оценить дисперсии погрешностей, то есть квадрат остатка -й выборки:

Необходимые остатки возвращаются при невзвешенной подгонке (OLS), так как оценка OLS является несмещенной. Затем передайте обратную величину остатков в качестве оценочного оптимального веса для выполнения метода взвешенных наименьших квадратов:

Hansen, P. C., Pereyra, V. and Scherer, G. (2013) Least squares data fitting with applications. Baltimore, Md: Johns Hopkins University Press, стр. 147—​155.

Kutner, M. H. et al. (2005) Applied Linear statistical models.

Weisberg, S. (2014) Applied linear regression. Fourth edition. Hoboken, NJ: Wiley (Wiley series in probability and statistics).