Метод аппроксимации в статистике
Нелинейный метод наименьших квадратов — это форма анализа наименьших квадратов , используемая для подгонки набора из m наблюдений с помощью модели, которая является нелинейной по n неизвестным параметрам ( m ≥ n ). Он используется в некоторых формах нелинейной регрессии . Основой метода является аппроксимация модели линейной и уточнение параметров последовательными итерациями. Существует много сходств с линейным методом наименьших квадратов , но также имеются и некоторые существенные различия . В экономической теории нелинейный метод наименьших квадратов применяется в (i) пробит-регрессии, (ii) пороговой регрессии, (iii) гладкой регрессии, (iv) логистической регрессии связей, (v) регрессорах, преобразованных по Боксу–Коксу ( ).
Теория
Рассмотрим набор точек данных и кривую (функцию модели) , которая в дополнение к переменной также зависит от параметров, причем требуется найти вектор параметров, такой, чтобы кривая наилучшим образом соответствовала заданным данным в смысле наименьших квадратов, то есть сумма квадратов
была бы минимизирована, где остатки (ошибки прогнозирования внутри выборки) r i определяются как
для
Минимальное значение S достигается, когда градиент равен нулю . Поскольку модель содержит n параметров, то существует n уравнений градиента:
В нелинейной системе производные являются функциями как независимой переменной, так и параметров, поэтому в общем случае эти градиентные уравнения не имеют замкнутого решения. Вместо этого для параметров должны быть выбраны начальные значения. Затем параметры уточняются итеративно, то есть значения получаются путем последовательного приближения,
Здесь k — номер итерации, а вектор приращений известен как вектор сдвига. На каждой итерации модель линеаризуется путем аппроксимации к разложению полинома Тейлора первого порядка вокруг
Матрица Якоби , J , является функцией констант, независимой переменной и параметров, поэтому она изменяется от одной итерации к другой. Таким образом, в терминах линеаризованной модели,
а остатки задаются как
Подставляя эти выражения в градиентные уравнения, они принимают вид
которые при перестановке становятся n одновременными линейными уравнениями, нормальными уравнениями
Нормальные уравнения записываются в матричной форме как
Эти уравнения составляют основу алгоритма Гаусса–Ньютона для нелинейной задачи наименьших квадратов.
Обратите внимание на соглашение о знаках в определении матрицы Якоби в терминах производных. Формулы, линейные по , могут появляться с множителем в других статьях или литературе.
Расширение с помощью весов
Если наблюдения не одинаково надежны, взвешенную сумму квадратов можно минимизировать,
Каждый элемент диагональной весовой матрицы W в идеале должен быть равен обратной величине дисперсии ошибки измерения. [1]
Нормальные уравнения тогда, в более общем виде, имеют вид:
Геометрическая интерпретация
В линейном наименьших квадратах целевая функция S является квадратичной функцией параметров.
Когда есть только один параметр, график S относительно этого параметра будет параболой . При двух или более параметрах контуры S относительно любой пары параметров будут концентрическими эллипсами (предполагая, что матрица нормальных уравнений положительно определена ) . Минимальные значения параметров находятся в центре эллипсов. Геометрию общей целевой функции можно описать как параболоид эллиптический. В NLLSQ целевая функция является квадратичной относительно параметров только в области, близкой к ее минимальному значению, где усеченный ряд Тейлора является хорошим приближением к модели.
Чем больше значения параметров отличаются от их оптимальных значений, тем больше контуры отклоняются от эллиптической формы. Следствием этого является то, что начальные оценки параметров должны быть как можно ближе к их (неизвестным!) оптимальным значениям. Также объясняется, как может возникнуть расхождение, поскольку алгоритм Гаусса–Ньютона сходится только тогда, когда целевая функция приблизительно квадратична по параметрам.
Вычисление
Первоначальные оценки параметров
Некоторые проблемы плохой обусловленности и расхождения можно исправить, найдя начальные оценки параметров, близкие к оптимальным значениям. Хороший способ сделать это — компьютерное моделирование . Как наблюдаемые, так и вычисленные данные отображаются на экране. Параметры модели корректируются вручную до тех пор, пока согласие между наблюдаемыми и вычисленными данными не станет достаточно хорошим. Хотя это будет субъективным суждением, достаточно найти хорошую отправную точку для нелинейного уточнения. Начальные оценки параметров можно создать с помощью преобразований или линеаризации. Еще лучше, эволюционные алгоритмы, такие как алгоритм стохастической воронки, могут привести к выпуклому бассейну притяжения, который окружает оптимальные оценки параметров. [ требуется цитата ] Гибридные алгоритмы, которые используют рандомизацию и элитизм, а затем методы Ньютона, как было показано, полезны и вычислительно эффективны [ необходима цитата ] .
Решение
Для поиска решения можно применить любой из описанных ниже методов.
Критерии сходимости
Критерий сходимости, основанный на здравом смысле, заключается в том, что сумма квадратов не увеличивается от одной итерации к другой. Однако этот критерий часто трудно реализовать на практике по разным причинам. Полезный критерий сходимости:
Значение 0,0001 является несколько произвольным и может нуждаться в изменении. В частности, его может потребоваться увеличить, когда экспериментальные ошибки велики. Альтернативный критерий:
Опять же, численное значение несколько произвольно; 0,001 эквивалентно указанию того, что каждый параметр должен быть уточнен до точности 0,1%. Это разумно, когда оно меньше наибольшего относительного стандартного отклонения параметров.
Расчет якобиана методом численного приближения
Существуют модели, для которых очень сложно или даже невозможно вывести аналитические выражения для элементов якобиана. Тогда численное приближение
получается путем вычисления для и . Размер приращения, , следует выбирать так, чтобы численная производная не подвергалась ошибке приближения, будучи слишком большой, или ошибке округления , будучи слишком маленькой.
Ошибки параметров, доверительные интервалы, остатки и т. д.
Некоторая информация приведена в соответствующем разделе на странице «Взвешенные наименьшие квадраты» .
Множественные минимумы
Множественные минимумы могут возникать при различных обстоятельствах, вот некоторые из них:
- Параметр возводится в степень двух или более. Например, при подгонке данных к кривой Лоренца , где — высота, — положение, а — полуширина на половине высоты, существуют два решения для полуширины, которые дают одно и то же оптимальное значение для целевой функции.
- Два параметра можно поменять местами, не меняя значения модели. Простой пример — когда модель содержит произведение двух параметров, так как даст то же значение, что и .
- Параметр находится в тригонометрической функции, например , , которая имеет одинаковые значения при . См. пример алгоритма Левенберга–Марквардта .
Не все множественные минимумы имеют одинаковые значения целевой функции. Ложные минимумы, также известные как локальные минимумы, возникают, когда значение целевой функции больше, чем ее значение в так называемом глобальном минимуме. Чтобы быть уверенным, что найденный минимум является глобальным минимумом, уточнение следует начинать с сильно различающимися начальными значениями параметров. Когда один и тот же минимум найден независимо от начальной точки, он, скорее всего, будет глобальным минимумом.
При наличии множественных минимумов возникает важное последствие: целевая функция будет иметь максимальное значение где-то между двумя минимумами. Матрица нормальных уравнений не является положительно определенной в максимуме целевой функции, так как градиент равен нулю и не существует уникального направления спуска. Уточнение из точки (набора значений параметров), близкой к максимуму, будет плохо обусловленным и его следует избегать в качестве отправной точки. Например, при подгонке лоренцева матрица нормальных уравнений не является положительно определенной, когда полуширина полосы равна нулю. [2]
Преобразование в линейную модель
Нелинейная модель иногда может быть преобразована в линейную. Такая аппроксимация, например, часто применима вблизи наилучшей оценки, и это одно из основных предположений в большинстве итеративных алгоритмов минимизации. Когда линейная аппроксимация действительна, модель может быть напрямую использована для вывода с обобщенным наименьшими квадратами , где применяются уравнения линейной подгонки шаблона [3] .
Другим примером линейной аппроксимации может быть случай, когда модель представляет собой простую экспоненциальную функцию,
которую можно преобразовать в линейную модель путем логарифмирования.
Графически это соответствует работе на полулогарифмическом графике . Сумма квадратов становится
Этой процедуры следует избегать, если только ошибки не являются мультипликативными и логарифмически нормально распределенными, поскольку она может дать вводящие в заблуждение результаты. Это происходит из-за того, что какими бы ни были экспериментальные ошибки по y , ошибки по log y различны. Поэтому, когда преобразованная сумма квадратов минимизируется, будут получены разные результаты как для значений параметров, так и для их вычисленных стандартных отклонений. Однако с мультипликативными ошибками, которые распределены логарифмически нормально, эта процедура дает несмещенные и последовательные оценки параметров.
Другим примером является кинетика Михаэлиса–Ментен , используемая для определения двух параметров и : График
Лайнуивера –Берка
относительно является линейным по параметрам и , но очень чувствителен к ошибкам в данных и сильно смещен в сторону подгонки данных в определенном диапазоне независимой переменной .
Алгоритмы
Метод Гаусса–Ньютона
Нормальные уравнения
могут быть решены с помощью разложения Холецкого , как описано в линейном методе наименьших квадратов . Параметры обновляются итеративно,
где k — номер итерации. Хотя этот метод может быть адекватным для простых моделей, он не сработает, если возникнет расхождение. Поэтому защита от расхождения имеет важное значение.
Сдвиг-резка
Если происходит расхождение, простым приемом является уменьшение длины вектора сдвига, , на дробь, f
Например, длина вектора сдвига может быть последовательно уменьшена вдвое до тех пор, пока новое значение целевой функции не станет меньше ее значения на последней итерации. Дробь, f может быть оптимизирована линейным поиском . [4] Поскольку каждое пробное значение f требует пересчета целевой функции, не стоит оптимизировать ее значение слишком строго.
При использовании shift-cutting направление вектора сдвига остается неизменным. Это ограничивает применимость метода ситуациями, когда направление вектора сдвига не сильно отличается от того, каким оно было бы, если бы целевая функция была приблизительно квадратичной по параметрам,
Параметр Марквардта
Если происходит расхождение и направление вектора сдвига настолько далеко от его «идеального» направления, что сокращение сдвига не очень эффективно, то есть доля, f, необходимая для избежания расхождения, очень мала, направление должно быть изменено. Этого можно достичь с помощью параметра Марквардта . [5] В этом методе нормальные уравнения модифицируются
, где — параметр Марквардта, а I — единичная матрица. Увеличение значения приводит к изменению как направления, так и длины вектора сдвига. Вектор сдвига поворачивается в направлении наискорейшего спуска , когда — вектор наискорейшего спуска. Таким образом, когда становится очень большим, вектор сдвига становится малой долей вектора наискорейшего спуска.
Для определения параметра Марквардта были предложены различные стратегии. Как и в случае со сдвигом-урезанием, слишком строго оптимизировать этот параметр расточительно. Вместо этого, как только найдено значение, которое приводит к уменьшению значения целевой функции, это значение параметра переносится на следующую итерацию, уменьшается, если это возможно, или увеличивается, если это необходимо. При уменьшении значения параметра Марквардта существует пороговое значение, ниже которого можно безопасно установить его равным нулю, то есть продолжить с немодифицированным методом Гаусса–Ньютона. Пороговое значение может быть установлено равным наименьшему сингулярному значению якобиана. [6] Граница для этого значения задается выражением , где tr — функция следа . [7]
QR-разложение
Минимум суммы квадратов можно найти методом, не включающим формирование нормальных уравнений. Остатки с линеаризованной моделью можно записать как
Якобиан подвергается ортогональному разложению; разложение QR будет служить для иллюстрации процесса.
где Q — ортогональная матрица, а R — матрица, которая разделена на блок, и нулевой блок. — верхний треугольный.
Остаточный вектор умножается слева на .
Это не влияет на сумму квадратов, поскольку Q ортогонален . Минимальное значение S достигается , когда верхний блок равен нулю. Поэтому вектор сдвига находится путем решения
Эти уравнения легко решаются, поскольку R — верхний треугольный элемент.
Разложение по сингулярным значениям
Вариант метода ортогонального разложения включает в себя разложение по сингулярным значениям , в котором R диагонализируется с помощью дальнейших ортогональных преобразований.
где ортогонален, — диагональная матрица сингулярных значений, а — ортогональная матрица собственных векторов или, что эквивалентно, правых сингулярных векторов . В этом случае вектор сдвига задается как
Относительная простота этого выражения очень полезна в теоретическом анализе нелинейных наименьших квадратов. Применение разложения сингулярных значений подробно обсуждается в Lawson и Hanson. [6]
Градиентные методы
В научной литературе имеется множество примеров использования различных методов для решения задач нелинейной подгонки данных.
- Включение вторых производных в разложение ряда Тейлора функции модели. Это метод Ньютона в оптимизации . Матрица H известна как матрица Гессе . Хотя эта модель имеет лучшие свойства сходимости вблизи минимума, она намного хуже, когда параметры далеки от своих оптимальных значений. Вычисление Гессе увеличивает сложность алгоритма. Этот метод не является общепринятым.
- Метод Дэвидона–Флетчера–Пауэлла . Этот метод, разновидность метода псевдоНьютона, похож на описанный выше, но вычисляет гессиан методом последовательного приближения, чтобы избежать необходимости использовать аналитические выражения для вторых производных.
- Наискорейший спуск . Хотя уменьшение суммы квадратов гарантируется, когда вектор сдвига указывает в направлении наискорейшего спуска, этот метод часто работает плохо. Когда значения параметров далеки от оптимальных, направление вектора наискорейшего спуска, которое является нормальным (перпендикулярным) к контурам целевой функции, сильно отличается от направления вектора Гаусса–Ньютона. Это делает расхождение гораздо более вероятным, особенно потому, что минимум вдоль направления наискорейшего спуска может соответствовать небольшой доле длины вектора наискорейшего спуска. Когда контуры целевой функции очень эксцентричны из-за высокой корреляции между параметрами, итерации наискорейшего спуска с отсечением сдвига следуют медленной зигзагообразной траектории к минимуму.
- Поиск сопряженного градиента . Это улучшенный метод наискорейшего спуска с хорошими теоретическими свойствами сходимости, хотя он может не работать на цифровых компьютерах с конечной точностью даже при использовании для квадратичных задач. [8]
Методы прямого поиска
Методы прямого поиска зависят от оценок целевой функции при различных значениях параметров и вообще не используют производные. Они предлагают альтернативы использованию числовых производных в методе Гаусса–Ньютона и градиентных методах.
- Поиск с чередующимися переменными. [4] Каждый параметр изменяется по очереди путем добавления к нему фиксированного или переменного приращения и сохранения значения, которое приводит к уменьшению суммы квадратов. Метод прост и эффективен, когда параметры не сильно коррелируют. Он имеет очень плохие свойства сходимости, но может быть полезен для поиска начальных оценок параметров.
- Поиск Нелдера–Мида (симплекс) . Симплекс в этом контексте — это многогранник из n + 1 вершин в n измерениях; треугольник на плоскости, тетраэдр в трехмерном пространстве и т. д. Каждая вершина соответствует значению целевой функции для определенного набора параметров. Форма и размер симплекса настраиваются путем изменения параметров таким образом, что значение целевой функции в наивысшей вершине всегда уменьшается. Хотя сумма квадратов может изначально быстро уменьшаться, она может сходиться к нестационарной точке на квазивыпуклых задачах, по примеру MJD Powell.
Более подробные описания этих и других методов доступны в книге «Численные рецепты » вместе с компьютерным кодом на разных языках.
Смотрите также
Ссылки
- ^ Это подразумевает, что наблюдения некоррелированы. Если наблюдения коррелированы , выражение применимо. В этом случае матрица весов в идеале должна быть равна обратной матрице дисперсии-ковариации ошибок наблюдений.
- ^ При отсутствии ошибки округления и экспериментальной ошибки в независимой переменной матрица нормальных уравнений была бы сингулярной
- ^ Britzger, Daniel (2022). "The Linear Template Fit". Eur. Phys. J. C. 82 ( 8): 731. arXiv : 2112.01548 . Bibcode : 2022EPJC...82..731B. doi : 10.1140/epjc/s10052-022-10581-w.
- ^ ab MJ Box, D. Davies и WH Swann, Нелинейные методы оптимизации, Oliver & Boyd, 1969
- ^ Этот метод был предложен независимо Левенбергом (1944), Жираром (1958), Уинном (1959), Моррисоном (1960) и Марквардтом (1963). В большинстве научных работ для него используется только имя Марквардта. Ссылки на цитирование см. в основной статье .
- ^ ab CL Lawson и RJ Hanson, Решение задач наименьших квадратов, Prentice–Hall, 1974
- ^ Р. Флетчер, UKAEA Report AERE-R 6799, HM Stationery Office, 1971
- ^ MJD Powell, Computer Journal, (1964), 7 , 155.
Дальнейшее чтение
- Келли, CT (1999). Итеративные методы оптимизации (PDF) . SIAM Frontiers in Applied Mathematics. Том № 18. ISBN 0-89871-433-8.
- Штрутц, Т. (2016). Подгонка данных и неопределенность: практическое введение в метод взвешенных наименьших квадратов и далее (2-е изд.). Springer Vieweg. ISBN 978-3-658-11455-8.