CC BY
116
НАУЧНОЕ ИЗДАНИЕ МГТУ ИМ. Н. Э. БАУМАНА
НАУКА и ОБРАЗОВАНИЕ
Эл № ФС77 - 48211. Государственная регистрация №0421200025. ISSN 1994-0408
электронный научно-технический журнал
Алгоритм вычисления М-оценок параметров авторегрессионного поля # 07, июль 2013 Б01:10.7463/0713.0571094 Горяинов В. Б.
УДК 519.234.3
Россия, МГТУ им. Н.Э. Баумана vb-goryainov@mail.ru
Введение
Рассмотрим модель пространственной авторегрессии — стационарное поле Ху, описываемое уравнением
Ху = атХг-1,у + аыХ^-! + ацХ*-^— + еу, г,] = 0, ±1, ±2,..., (1)
где а = (аю, а01, ац)т — авторегрессионные коэффициенты, а еу — независимые одинаково распределенные случайные величины с нулевым математическим ожиданием Ееу = 0. Такие поля широко используются в естественных, технических и гуманитарных науках (см., например, [1,2,3,4]).
Одной из важных задач является построение оценок коэффициентов а по наблюдениям авторегрессионного поля Ху. Наиболее распространенным методом оценивания до сих пор является метод максимального правдоподобия [5]. Метод максимального правдоподобия предполагает известным тип функции распределения ^(ж) обновляющего поля еу. В частности, если еу имеют нормальное распределение, то оценки максимального правдоподобия совпадают с оценками наименьших квадратов и с асимптотически эквивалентными им оценками типа Юла — Уолкера [6]. Если предположения о типе ^(ж) выполняются абсолютно точно, то метод максимального правдоподобия является оптимальным. Однако при решении практических задач тип распределения ^(ж) обновляющего поля еу известен лишь приближенно или неизвестен совсем. В этом случае оценки максимального правдоподобия теряют свою эффективность [7].
Одной из альтернатив оценкам максимального правдоподобия являются М-оценки (оценки, обобщающие оценки максимального правдоподобия). М-оценки хорошо зарекомендовали себя в более простых моделях (одно- и двухвыборочная задача и линейная регрессия [8,9],
процессы авторегрессии-скользящего среднего [10]). В этих моделях М-оценки оказались робастными, т.е. слабо чувствительными к нарушению предположения о виде распределения наблюдений. Свойства М-оценок коэффициентов авторегрессионного уравнения (1) изучались в [7], где была доказана асимптотическая нормальность М-оценок и вычислена асимптотическая относительная эффективность М-оценок по отношению к оценкам наименьших квадратов. Там же было показано, что уже при небольшом нарушении предположения о нормальности е^ использование М-оценок предпочтительнее оценок наименьших квадратов.
Нахождение М-оценки коэффициентов а в модели (1) сводится к решению системы уравнений, для чего формально можно использовать метод Ньютона (также известный как метод касательных) — итерационный численный метод нахождения нуля заданной функции. Однако специфика этой системы уравнений делает процедуру вычисления М-оценок методом метод касательных Ньютона ненадежной. Возникает потребность в алгоритме, который является заведомо сходящимся.
В данной работе построен алгоритм вычисления М-оценок коэффициентов а в модели (1) и доказана его сходимость. Алгоритм представляет собой итерационный вариант взвешенного метода наименьших квадратов с пересчитываемыми на каждом шаге весами и основан на идее из [11].
1. Определение М-оценок
Рассмотрим поле (1), где а = (аю, а01, ац) — неизвестный вектор параметров, а е^ — независимые одинаково распределенные случайные величины с плотностью f (ж). Как показано в [6], достаточным условием стационарности Х^ является отсутствие корней уравнения
1 — аю^1 — ао1^2 — = 0
внутри единичного полидиска |г1| < 1, \г2\ < 1, что равносильно выполнению условий [12]
\аря1 < 1 (РЛ) еХ;
/ О О 0\0 / \0
(1 + а10 — аС)1 — а11) — 4(аш + ао1аи) > 0;
1 — а01 > \аю + ао1ап\.
ЗдесьХ = {(0,1), (1,0), (1,1)}.
М-оценку а*тп параметра а по наблюдениям Х^, г = 1,..., т, ] = 1,..., п, поля (1) определим как точку минимума функции
т п
^(а) = ^^ ^^ Р(ХгЗ — a10Xi—1,j — a01Xi,j-1 — a11Xi-1,j-1), (2)
г=2 j=2
где р — подходящим образом выбранная функция. Обычно р(ж) — выпуклая четная гладкая функция.
В [7] доказана следующая теорема.
Теорема 1. Пусть поле Х^, описываемое уравнением (1), является стационарным, р(ж) — выпуклая функция, р''(ж) непрерывна почти всюду и ограничена, а плотность / (ж) независимых одинаково распределенных случайных величин е^ в (1) удовлетворяет следующим условиям: Беи = 0, Ер'(ец) = 0, Е(р'(ец)2) > 0, а2 = йен < то, 0 < Ер"(ец) < то.
Тогда при т, п ^ то случайный вектор \]тп(атп — а) является асимптотически нормальным с нулевым математическим ожиданием и ковариационной матрицей
,-х Е[р'(ец)2]
Я_
(Е[р"(еп)])2:
где Я — ковариационная матрица вектора (Хю, Х01, Хц )т.
Заметим, что если / (ж) — четная функция (а это типичная ситуация), то производная р'(ж) четной функции р(ж) является нечетной и условие Ер'(е11) = 0 будет выполнено.
В [13] асимптотическая нормальность М-оценок и их обобщений доказана и для невыпуклой функции р(ж), если только существует р/;/(ж).
Если функция р(ж) является выпуклой и дифференцируемой, то минимизация (2) равносильна решению системы уравнений
Лрд(а) = 0, р,д ЕХ, (3)
где
т п
Лрд (а) = ^ ^ ф(Х^ — аюХ^ц — ао1Х^_1 — аиХ^ц_1)Х^ра_д = 0, р, д Е Х, (4)
г=2 j=2
а ф(ж) = р'(ж).
Из (1) следует, что функция правдоподобия параметра а имеет вид
тп
^(а) = ^ /(Xij — а10Xi_ 1,j — a01Xi,j_1 — a11Xi_1,j_1). i=2 j=2
Оценка максимального правдоподобия параметра а определяется как точка максимума Ь(а), или, что равносильно, как точка минимума функции
тп
1(а) = — 1п £(а) = ^^ (— 1п / (Х^ — аюХ^^- — ао1Х^_1 — ацХ^,^)). i=2 j=2
Таким образом, оценка максимального правдоподобия является частным случаем М-оценки с р-функцией р(ж) = — 1п /(ж), или ф-функцией
/'(ж) / (ж)
Рекомендации по выбору р(ж) имеются в [8,9]. Наиболее распространенной р-функцией является р-функция Хьюбера — параметрическое семейство, задаваемое уравнением
Мж2 I жI <СГ к*
\ I ? | I ' ' 7
ж) = <
; 2к|ж| — к2, |ж| >к.
Из определения р-функции Хьюбера следует, что если невязка
— а10Хг-1,.? — а01Хг^-1 — а11Х1-1,]-1 (5)
в минимизируемой сумме (2) принимает достаточно большие значения, то ее вклад в эту сумму линеен, т.е. ограничивается по сравнению с методом наименьших квадратов. Если же эта величина принимает небольшие значения, то ее вклад в сумму (2) такой же, как и в методе наименьших квадратов.
Еще в большей степени минимизировать влияние больших значений невязки (5) на сумму (2) можно, взяв в качестве рТ(х) функцию, называемую бивес Тьюки:
РТ(х) = { 1 - 0 - (Х)2|х| £ к
I 1, |х| > к.
Из определения р-функции Тьюки видно, что если значение невязки (5) становится больше, чем значение параметра к, то ее вклад в сумму (2) практически игнорируется, заменяясь единицей.
Таким образом, по сравнению с оценками наименьших квадратов, М-оценки с р-функция-ми Хьюбера и Тьюки должны быть менее чувствительными к большим значениям слагаемых в (2), которые индуцируются большими ошибками в наблюдениях поля Х^. Смысл понятий «небольшие значения» и «достаточно большие значения» регулируется параметром к. В частности, если р(х) = х2, то получается оценка наименьших квадратов, а если р(х) = |х| — оценка наименьших модулей.
2. Построение алгоритма
При решении системы уравнений (3) формально можно использовать алгоритм Ньютона [14], суть которого состоит в линеаризации уравнений этой системы, т.е. в замене функций ЛРд (а) их разложением по формуле Тейлора в точке текущего приближения решения. Таким образом, если к-ю аппроксимацию решения обозначить а*к, то (к + 1)-я итерация является решением системы уравнений
л/ I дЛ(ак)(* *\ п
Л(ак) + да (ак+1 - ак) =
К сожалению, этот алгоритм может расходиться, например, при неудачном выборе начально-
« д Л(а*к)
го приближения. Кроме того, если определитель матрицы — к может быть сколь угодно
малым, то скорость сходимости последовательности а*к не будет квадратичной, а сам метод может преждевременно прекратить поиск и дать неверное для заданной точности приближение.
Для того чтобы М-оценки были робастными (устойчивыми к нарушению предположения о распределении обновляющего поля е^), функция ^(х) в (4) должна быть ограничена.
В этом случае ее производная ^'(х) будет стремиться к нулю на бесконечности, и, следо-
дЛ(ак) _
вательно, определитель матрицы — к не будет отделен от нуля, что делает процедуру
вычисления М-оцеиок методом Ньютона ненадежной. По этой причине построим алгоритм, который являются заведомо сходящимся.
Предположим, что функция ф дифференцируема. Определим функцию
( ^, х = 0; т(х) = { ж
[ ф'(0), х = 0.
В этом случае решение системы (3) равносильно решению системы
т п
^^Цеу (а)) ец(а) Ха/3 = 0, а, в е!, (6)
г=2 3=2
где е3з (а) = Ху — а10Хг_133 — а01Хг>3_1 — ацХг-1з-1. Запишем систему (6) в векторном виде:
тп
ш(е»3 (а)) егз (а) Х^ = 0, (7)
i=2 3=2
где Хгз = (Хг_1,3-, Хг,3_1, Хг_1,3_1). Если бы величины и>г3- = и)(егз(а)) были известны (не зависели от а), то на решение системы (7) относительно а можно было бы смотреть как на оценки взвешенного метода наименьших квадратов, поскольку они получались бы как точка минимума функции
тп
/ у / у шг3 ег3(
г=2 3=2
(а) = ^ У" 1г3 е3(а). (8)
Отметим, что можно было бы свести минимизацию (8) и к нахождению обычных оценок наименьших квадратов, если в качестве «наблюдений» взять /1г3Хг3-, а в качестве векторов-столбцов регрессионной матрицы величины //Шг3"Хг3-. Тогда решение (7) — точка минимума функции
тп
Ьт (а) = ^^ Хг3 — Щ^ 2.
г=2 3=2
Эти соображения наводят на мысль нахождения минимума (2) итерационным способом. Соответствующий алгоритм будет иметь следующий вид.
1. Определяется начальное приближение а** М-оценки а* параметра а.
2. По известному к-му приближению а*к вычисляются
ег3к = Хг3 — а*кТ Хг3, 1г3к = 1(ег3к), г = 2,... ,т, ] = 2 ,...,п.
3. Вычисляется следующее (к + 1)-е приближение ак+1 как решение системы линейных относительно а уравнений
тп
(Хг3 — аТ Хг3 )Хг3 = 0
г=2 3=2
или
т п т п
£ £ Хг3 ХТа = £ £ 13 (Хг3Хг3). (9)
г=2 3=2 г=2 3=2
4. Процесс заканчивается, когда
К+1- а*1 < 6
где постоянная 6 задает точность вычисления М-оценки а*. Обозначив
Л,(а) = ^ Це^(а)) (Х^ - вТХу)
в г=2 ]=2
алгоритм можно записать в виде
ак+1 = Мак), (10)
причем решение а системы (3) удовлетворяет условию
=
3. Сходимость алгоритма
Докажем сходимость этого алгоритма.
Теорема 2. Пусть функция р(х) не убывает и имеет производную р'(х), а эд(х) не возрастает при |х| ^ то. Пусть существует единственное решение а * системы (3). Тогда алгоритм сходится к а .
Доказательство. Определим функцию д(х) по формуле
д(х) = р(^х).
Тогда р(х) = д(х2), р'(х) = 2хд'(х2) и эд(х) = 2д'(ж2) при х = 0. Применив теорему о среднем, получим
т п т п
ь(ак+1) - вд) = XX (д(4 (ак+1)) - д(4 (ак )Й = XX д'(т К4 (ак+1) - 4 (ак Й,
г=2 .7=2 г=2 .7=2
где т лежит между (ак) и (а£+1).
Заметим, что если и>(х) не возрастает при х > 0, то д'(х) также не возрастает. Поэтому д'(т) < д' (х) при т > х. Следовательно, если х < у, то
дЫ - д(х) = д'(т)(у - х) < д'(х)(у - x), х < т < у.
Если же х > у, то
д(х) - д(у) = д'(т)(х - у) > д'(х)(х - у), х < т < у
и поэтому также
д(у) - д(х) = д'(т)(у - х) < д'(х)(у - x), х < т < у.
Таким образом,
тп
ь(ак+1) - ь(а1) < X X д' (4 (ак )К 4 (ак+1) - 4 (ак Й.
г=2 ¿=2
Так как
4 (4+1) — 4 (4) = (4 — 4+0 "г.' 4 (4+1) + 4 (4) = — (4 + 4+1) "г. ,
то
¿(4+1) - ¿(аЮ < Е ^'(4 (ак й (4+1)- % (4 й (ак+1)+(4))
г=2 .=2
= Е Е я'(4 (акЙ К - 4+1)Т (2Ху - (4 + 4+1 .
г=2 .=2
Так как
0'(4 К)) = 1 ™(4 (ак)),
1
и ак+1 является решением (9), то
ЕЕи(4 (ак Й аТ+1 -"г. = ЕЕ и(4 (ак Й .
г=2 .=2 г=2 .=2
Поэтому
1 т п
¿(4+1) - ь(ак) < ^ Е Е(ак- ак+1)Т -"г. 44 (ак й (2ак+А- - (ак + ак+1)Т ".О =
¿=2 .=2
1 т п
=1 ЕЕ и(4 (ак Й(4- ак+1)Т "г. (4+1- 4).
2 г=2 .=2
Матрица X. X. положительно определена. Следовательно,
(4 — ак+1)Т"г.""г. (ак+1 — 4) < 0.
Так как р(х) возрастает при ж > 0, то ^(х) > 0 при ж > 0 и, следовательно, и(х) > 0 при х > 0. Поэтому (а£)) > 0 и, следовательно, Ь(а^+1) — ¿(ак) < 0.
Докажем сходимость последовательности а£ к М-оценке а* параметра а. Так как последовательность ¿(ак) не возрастает и ограничена снизу, то она имеет конечный предел. Следовательно, последовательность а£ ограничена. Действительно, в противном случае существовала бы подпоследовательность , стремящаяся к бесконечности, и тогда подпоследовательность ) в силу возрастания функции р(х) при х > 0 также стремилась бы к бесконечности.
Так как а£ ограничена, то у нее существует сходящаяся подпоследовательность, предел которой обозначим а^. Перейдя в обеих частях равенства (10) к пределу, получим, что а^> = Л(а^), т.е. а^ — стационарная точка преобразования Л.
Если а^ — единственная стационарная точка Л, то а к ^ а^. Действительно, в противном случае существовала бы подпоследовательность 6^. последовательности а£, для которой для некоторого 8 > 0 |6£. — а^| > 8. Последовательность 6^. в свою очередь имела бы сходящуюся подпоследовательность, стремящуюся к 6^ = а^, и для которой бы выполнялось бы 6^ = Л(6^). Теорема доказана.
4. Выводы
В работе построен алгоритм вычисления М-оценок коэффициентов уравнения авторегрессионного поля и доказана сходимость этого алгоритма. Алгоритм представляет собой итерационный вариант взвешенного метода наименьших квадратов с пересчитываемыми на каждом шаге весами. В отличие от метода Ньютона (метода касательных) алгоритм является сходящимся при любом выборе начального приближения.
Список литературы
1. Ripley B.D. Spatial statistics. Hoboken: Wiley, 1981. DOI: 10.1002/0471725218.
2. Gaetan C., Guyon X. Spatial statistics and modeling. New York: Springer, 2010. (Springer Series in Statistics.) DOI: 10.1007/978-0-387-92257-7.
3. Handbook of spatial statistics / A.E. Gelfand, P.J. Diggle, M. Fuentes, P. Guttorp (eds.). Boca Raton: Taylor & Francis, 2010.
4. HainingR. Spatial Data Analysis: Theory and Practice. Cambridge: CUP, 2004.
5. Yao Q., Brockwell P.J. Gaussian Maximum Likelihood Estimation for ARMA Models II Spatial Processes // Bernoulli. 2006. Vol. 12, no. 3. P. 403-429.
6. Tjostheim D. Statistical Spatial Series Modelling // Advances in Applied Probability. 1978. Vol. 10, no. 1.P. 130-154.
7. Горяинов В.Б. М-оценки пространственной авторегрессии // Автоматика и телемеханика. 2012. №8. С. 119-129.
8. Хьюбер П.Дж. Робастность в статистике. М.: Мир, 1984.
9. Хампель Ф., Рончетти Э., Рауссеу П., Штаэль В. Робастность в статистике. Подход на основе функций влияния. М.: Мир, 1989.
10. Lee C.-H., Martin R.D. Ordinary and proper location M-estimates for autoregressive-moving average models //Biometrika. 1986. Vol. 73, no. 3. P. 679-686.
11. Stirling W.D. Iteratively reweighted least squares for models with a linear part // J. Roy. Statist. Soc. Ser. C. 1984. Vol. 33, no. 1. P. 7-17.
12. Basu S., Reinsel G.C. Properties ofthe Spatial Unilateral First-Order ARMA Model //Advances in Applied Probability. 1993. Vol. 25, no. 3. P. 631-648.
13. Горяинов В.Б. Обобщенные М-оценки коэффициентов авторегрессионного поля .. Автоматика и телемеханика. 2012. № 10. С. 42-51.
14. Амосов А.А., Дубинский Ю.А., Копченова Н.В. Вычислительные методы для инженеров. М.: Высшая школа, 1994.
SCIENTIFIC PERIODICAL OF THE BAUMAN MSTU
SCIENCE and EDUCATION
EL № FS77 - 48211. №0421200025. ISSN 1994-0408
electronic scientific and technical journal
Algorithm for calculating M-estimates of autoregressive field's parameters # 07, July 2013 DOI: 10.7463/0713.0571094 Goryainov V. B.
Bauman Moscow State Technical University 105005, Moscow, Russian Federation vb-goryainov@mail.ru
A process of a two-dimensional autoregression of order (1, 1) is considered in this article. Distribution of the innovation field of the autoregressive model is assumed to be unknown. An algorithm for calculating M-estimates of coefficients of the autoregressive field was constructed. The convergence of this algorithm was proved. The algorithm is an iterative version of the weighted least squares method. The weights are recalculated at each step. Each iteration represents the process of solving the system of linear equations. In contrast to the method of Newton (method of tangents) the algorithm converges from any point of the initial approximation.
References
1. Ripley B.D. Spatial statistics. Hoboken, Wiley, 1981. DOI: 10.1002/0471725218.
2. Gaétan C., Guyon X. Spatial statistics and modeling. Springer New York, 2010. (Springer Series in Statistics). DOI: 10.1007/978-0-387-92257-7.
3. Gelfand A.E., Diggle P.J., Fuentes M., Guttorp P. (eds.). Handbook of spatial statistics. Boca Raton: Taylor & Francis, 2010.
4. HainingR. Spatial Data Analysis: Theory and Practice. Cambridge, CUP, 2004.
5. Yao Q., Brockwell P.J. Gaussian Maximum Likelihood Estimation for ARMA Models II Spatial Processes. Bernoulli, 2006, vol. 12, no. 3, pp. 403-429.
6. Tjostheim D. Statistical Spatial Series Modelling. Advances in Applied Probability, 1978, vol. 10, no. 1, pp. 130-154.
7. Goriainov V.B. M-otsenki prostranstvennoi avtoregressii [M-estimates of the spatial autoregression coefficients]. Avtomatika i telemekhanika, 2012, no. 8, pp. 119-129. (Trans.
version: Automation and Remote Control, 2012, vol. 73, no. 8, pp. 1371-1379. DOI: 10.1134/S0005117912080103.)
8. HuberP.J. Robust Statistics. John Wiley & Sons, 1981. 320p. (Russ. ed.: HuberP.J. Robastnost v statistice. Moscow, Mir, 1984. 304 p.).
9. Hampel F.R., Ronchetti E.M., Rausseu P.J., Stahel W.A. Robust statistics: The approach based on influence functions. Wiley, New York, 1986. 536 p. (Wiley Series in Probability and Statistics). (Russ. ed.: Khampel' F., Ronchetti E., Rausseu P., Shtael' V. Robastnost vstatistike. Podhodna osnove funkcii vlujanija. Moscow, Mir, 1989. 512 p.).
10. Lee C.-H., Martin R.D. Ordinary and proper location M-estimates for autoregressive-moving average models. Biometrika, 1986, vol. 73, no. 3, pp. 679-686.
11. Stirling W.D. Iteratively reweighted least squares for models with a linear part. J. Roy. Statist. Soc. Ser. C, 1984, vol. 33, no. 1, pp. 7-17.
12. Basu S., Reinsel G.C. Properties of the Spatial Unilateral First-Order ARMA Model. Advances in Applied Probability, 1993, vol. 25, no. 3, pp. 631-648.
13. Goriainov V.B. Obobshchennye M-otsenki koeffitsientov avtoregressionnogo polia [Generalized M-estimates of the autoregression field coefficients]. Avtomatika i telemekhanika, 2012, no. 10, pp. 42-51. (Trans. version: Automation and Remote Control, 2012, vol. 73, no. 10, pp. 1624-1631. DOI: 10.1134/S0005117912100049.)
14. Amosov A.A., Dubinskii Iu.A., Kopchenova N.V. Vychislitel'nye metody dlia inzhenerov [Computational Methods for Engineers]. Moscow, Vysshaia shkola, 1994. 544 p.
