Ковариационные матрицы ошибок и оценок параметров эконометрических моделей

В предыдущих параграфах данной главы мы рассматривали случай отказа от одной предпосылки классической линейной модели множественной регрессии. Теперь мы расширим наш анализ и откажемся сразу от двух предпосылок: от предпосылок №№4-5 (о постоянстве дисперсии случайной ошибки и о некоррелированности разных случайных ошибок между собой). В техническом смысле этот параграф несколько сложнее предыдущих (в частности, тут более широко используется линейная алгебра). Поэтому, если вы заинтересованы в том, чтобы разобраться только в прикладных аспектах множественной регрессии, а в соответствующих вычислениях готовы полностью довериться эконометрическому пакету, можете его пропустить.

Отказ от указанных двух предпосылок означает, что ковариационная матрица вектора случайных ошибок (таблица, в которой записаны все ковариации между (varepsilon_{i}) и (varepsilon_{j}), см. параграф 3.3) больше не является диагональной матрицей с одинаковыми числами на главной диагонали, как это было в первоначальной классической модели. Теперь ковариационная матрица вектора случайных ошибок Ω — это произвольная ковариационная матрица (разумеется, так как это не совсем любая матрица, а именно ковариационная матрица, то по своим свойствам она является симметричной и положительно определенной).

Модель, в которой сохранены только первые три предпосылки классической линейной модели множественной регрессии, называется обобщенной линейной моделью множественной регрессии.

Проанализируем, к каким последствиям приводит отказ от предпосылок №№4-5.

Во-первых, полученная обычным методом наименьших квадратов оценка ({widehat{beta} = left( {X^{‘}X} right)^{- 1}}X’y) остается несмещенной (это свойство мы доказывали, опираясь как раз лишь на первые три предпосылки).

Во-вторых, МНК-оценки хоть и остаются несмещенными, но больше не являются эффективными.

В-третьих, если мы оцениваем ковариационную матрицу вектора оценок коэффициентов (которая нужна для тестирования всевозможных гипотез), то оценка (widehat{V}{{(widehat{beta})} = left( {X^{‘}X} right)^{- 1}}S^{2}) смещена и больше не является корректной.

Чтобы убедиться в этом, посчитаем ковариационную матрицу от (widehat{beta}) в условиях обобщенной модели (при этом мы используем свойства ковариационной матрицы, перечисленные в параграфе 3.3):

(V{left( widehat{beta} right) = V}{leftlbrack {{({X^{‘}X})}^{- 1}X’y} rightrbrack = V}{leftlbrack {{({X^{‘}X})}^{- 1}X'{({mathit{Xbeta} + varepsilon})}} rightrbrack =})

({}V{leftlbrack {{({{({X^{‘}X})}^{- 1}X’})}varepsilon} rightrbrack = left( {X^{‘}X} right)^{- 1}}X^{‘}Vlbrackvarepsilonrbrack{left( {left( {X^{‘}X} right)^{- 1}X^{‘}} right)^{‘} = {({X^{‘}X})}^{- 1}}X^{‘}Omega{left( {left( {X^{‘}X} right)^{- 1}X^{‘}} right)^{‘} = {({X^{‘}X})}^{- 1}}X^{‘}Omega X{({X^{‘}X})}^{- 1})

Так выглядит ковариационная матрица вектора МНК-оценок в обобщенной модели. Ясно, что она не может быть корректно оценена стандартной оценкой (left( {X^{‘}X} right)^{- 1}S^{2}). Следовательно, прежней формулой пользоваться нельзя: если мы будем использовать стандартные ошибки, рассчитанные по обычной формуле (предполагая выполнение предпосылок классической линейной модели), то получим некорректные стандартные ошибки, что может привести нас к неверным выводам по поводу значимости или незначимости тех или иных регрессоров.

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

Поэтому для получения эффективных оценок обычный МНК не подойдет, и придется воспользоваться альтернативным методом — обобщенным МНК (ОМНК, generalized least squares, GLS). Формулу для расчета оценок коэффициентов при помощи ОМНК позволяет получить специальная теорема.

Теорема Айткена

Если

1. модель линейна по параметрам и правильно специфицирована
   ({y = {mathit{Xbeta} + varepsilon}},)
2. матрица Х — детерминированная матрица, имеющая максимальный ранг k,
3. (E{(varepsilon) = overrightarrow{0}}),
4. (V{(varepsilon) = Omega}) — произвольная положительно определенная и симметричная матрица,

то оценка вектора коэффициентов модели ({{widehat{beta}}^{} = {({X’Omega^{- 1}X})}^{- 1}}X’Omega^{- 1}y) является:

• несмещенной
• и эффективной, то есть имеет наименьшую ковариационную матрицу в классе всех несмещенных и линейных по y оценок.

Предпосылки теоремы Айткена — это предпосылки обобщенной линейной модели множественной регрессии. Из них первые три — стандартные, как в классической модели, а четвертая ничего особого не требует (у вектора случайных ошибок может быть любая ковариационная матрица без каких-либо дополнительных специальных ограничений). Сама теорема Айткена является аналогом теоремы Гаусса — Маркова для случая обобщенной модели.

Докажем эту теорему.

Из линейной алгебры известно: если матрица (Omega) симметрична и положительно определена, то существует такая матрица P, что

(P’bullet{P = Omega^{- 1}}) ⇔ ({P^{‘} = Omega^{- 1}}bullet P^{- 1})

А раз такое представление возможно, то воспользуемся им для замены переменных. От вектора значений зависимой переменной (y), перейдем к вектору (Pbullet y), обозначив его как вектор ({y^{} = P}bullet y). Аналогичным образом введем матрицу ({X^{} = P}bullet X) и вектор ошибок ({varepsilon^{} = P}bulletvarepsilon).

Вернемся к исходной модели, параметры которой нас и интересуют:

({y = X}{beta + varepsilon})

Умножим левую и правую части равенства на матрицу (P):

({mathit{Py} = mathit{PX}}{beta + mathit{Pvarepsilon}})

С учетом новых обозначений это равенство можно записать так:

({y^{} = X^{}}{beta + varepsilon^{}})

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

Для этого вычислим математическое ожидание нового вектора ошибок:

(E{left( varepsilon^{} right) = E}{left( mathit{Pvarepsilon} right) = P}bullet E{(varepsilon) = P}bullet{overrightarrow{0} = overrightarrow{0}})

Теперь вычислим ковариационную матрицу вектора (varepsilon^{}):

(V{left( varepsilon^{} right) = V}{left( {Pbulletvarepsilon} right) = P}bullet V(varepsilon)bullet{P^{‘} = P}bulletOmegabullet{P^{‘} = P}bulletOmegabulletOmega^{- 1}bullet{P^{- 1} = I_{n}})

Здесь (I_{n}) обозначает единичную матрицу размера n на n.

Следовательно, для модели со звездочками выполняются все предпосылки теоремы Гаусса — Маркова. Поэтому получить несмещенную и эффективную оценку вектора коэффициентов можно, применив к этой измененной модели обычный МНК:

({{widehat{beta}}^{} = left( {{X^{}}^{‘}X^{}} right)^{- 1}}{X^{}}^{‘}y^{})

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

({{widehat{beta}}^{} = left( {{X^{}}^{‘}X^{}} right)^{- 1}}{X^{}}^{‘}{y^{} = left( {{(mathit{PX})}^{‘}mathit{PX}} right)^{- 1}}left( mathit{PX} right)^{‘}{mathit{Py} = left( {X’P’mathit{PX}} right)^{- 1}}X^{‘}P^{‘}{mathit{Py} =}left( {X’Omega^{- 1}P^{- 1}mathit{PX}} right)^{- 1}X’Omega^{- 1}P^{- 1}{mathit{Py} = {({X’Omega^{- 1}X})}^{- 1}}X’Omega^{- 1}y)

Что и требовалось доказать.

Взвешенный МНК, который мы обсуждали ранее, — это частный вариант обобщенного МНК (для случая, когда только предпосылка №4 нарушена, а предпосылка №5 сохраняется).

Как и при использовании взвешенного МНК в ситуации применения ОМНК коэффициент R-квадрат не обязан лежать между нулем и единицей и не может быть интерпретирован стандартным образом.

Слабая сторона ОМНК состоит в том, что для его реализации нужно знать не только матрицу регрессоров X с вектором значений зависимой переменной y, но и ковариационную матрицу вектора случайных ошибок (Omega). На практике, однако, эта матрица почти никогда не известна. Поэтому в прикладных исследованиях практически всегда вместо ОМНК используется, так называемый, доступный ОМНК (его ещё называют практически реализуемый ОМНК, feasible GLS). Идея доступного ОМНК состоит в том, что следует сначала оценить матрицу (Omega) (традиционно обозначим её оценку (widehat{Omega})), а уже затем получить оценку вектора коэффициентов модели, заменив в формуле ОМНК (Omega) на (widehat{Omega}):

({{widehat{beta}}^{} = {({X'{widehat{Omega}}^{- 1}X})}^{- 1}}X'{widehat{Omega}}^{- 1}y.)

Применение этого подхода осложняется тем, что (widehat{Omega}) не может быть оценена непосредственно без дополнительных предпосылок, так как в ней слишком много неизвестных элементов. Действительно, в матрице размер (n) на (n) всего (n^{2}) элементов, и оценить их все, имея всего (n) наблюдений, представляется слишком амбициозной задачей. Даже если воспользоваться тем, что матрица (Omega) является симметричной, в результате чего достаточно оценить только элементы на главной диагонали и над ней, мы все равно столкнемся с необходимостью оценивать (left( {n + 1} right){n/2}) элементов, что всегда больше числа доступных нам наблюдений.

Поэтому процедура доступного ОМНК устроена так:

1. Делаются некоторые предпосылки по поводу того, как устроена ковариационная матрица вектора случайных ошибок (Omega). На основе этих предпосылок оценивается матрица (widehat{Omega}).
2. После этого по формуле ({({X'{widehat{Omega}}^{- 1}X})}^{- 1}X'{widehat{Omega}}^{- 1}y) вычисляется вектор оценок коэффициентов модели.

Из сказанного следует, что доступный ОМНК может быть реализован только в ситуации, когда есть разумные основания сформулировать те или иные предпосылки по поводу матрицы (widehat{Omega}). Рассмотрим некоторые примеры таких ситуаций.

Пример 5.4. Автокорреляция и ОМНК-оценка

Рассмотрим линейную модель ({y = X}{beta + varepsilon}), для которой дисперсия случайных ошибок постоянна, однако наблюдается так называемая автокорреляция первого порядка:

(varepsilon_{i} = {{rho ast varepsilon_{i — 1}} + u_{i}})

Здесь (u_{i}) — независимые и одинаково распределенные случайные величины с дисперсией (sigma_{u}^{2}), а (rhoin{({{- 1},1})}) — коэффициент автокорреляции.

(а) Найдите ковариационную матрицу вектора случайных ошибок для представленной модели.

(б) Запишите в явном виде формулу ОМНК-оценки вектора коэффициентов модели, предполагая, что коэффициент (rho) известен.

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

Решение:

(а) Используя условие о постоянстве дисперсии случайной ошибки, то есть условие (mathit{var}{{(varepsilon_{i})} = mathit{var}}{(varepsilon_{i — 1})}), найдем эту дисперсию:

(mathit{var}{{(varepsilon_{i})} = mathit{var}}{({{rho ast varepsilon_{i — 1}} + u_{i}})})

(mathit{var}{{(varepsilon_{i})} = rho^{2}}mathit{var}{{(varepsilon_{i — 1})} + mathit{var}}{(u_{i})})

(mathit{var}{{(varepsilon_{i})} = rho^{2}}mathit{var}{{(varepsilon_{i})} + sigma_{u}^{2}})

(mathit{var}{left( varepsilon_{i} right) = frac{sigma_{u}^{2}}{1 — rho^{2}}}.)

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

(mathit{cov}{left( {varepsilon_{i},varepsilon_{i — 1}} right) = mathit{cov}}{left( {{{rho ast varepsilon_{i — 1}} + u_{i}},varepsilon_{i — 1}} right) = {rho ast mathit{cov}}}{left( {varepsilon_{i — 1},varepsilon_{i — 1}} right) + mathit{cov}}{left( {u_{i},varepsilon_{i — 1}} right) = {rho ast mathit{var}}}{{left( varepsilon_{i — 1} right) + 0} = frac{rho ast sigma_{u}^{2}}{1 — rho^{2}}})

По аналогии легко убедиться, что

(mathit{cov}{left( {varepsilon_{i},varepsilon_{i — k}} right) = frac{rho^{k} ast sigma_{u}^{2}}{1 — rho^{2}} = frac{sigma_{u}^{2} ast rho^{k}}{1 — rho^{2}}}.)

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

({Omega = frac{sigma_{u}^{2}}{1 — rho^{2}}}begin{pmatrix} begin{matrix} 1 & rho \ rho & 1 \ end{matrix} & ldots & begin{matrix} rho^{n — 2} & rho^{n — 1} \ rho^{n — 3} & rho^{n — 2} \ end{matrix} \ ldots & ldots & ldots \ begin{matrix} rho^{n — 2} & rho^{n — 3} \ rho^{n — 1} & rho^{n — 2} \ end{matrix} & ldots & begin{matrix} 1 & rho \ rho & 1 \ end{matrix} \ end{pmatrix})

(б) Вектор ОМНК-оценок коэффициентов имеет вид:

({{widehat{beta}}^{} = {({X’Omega^{- 1}X})}^{- 1}}X^{‘}Omega^{- 1}{y =}{left( {X’begin{pmatrix} begin{matrix} 1 & rho \ rho & 1 \ end{matrix} & ldots & begin{matrix} rho^{n — 2} & rho^{n — 1} \ rho^{n — 3} & rho^{n — 2} \ end{matrix} \ ldots & ldots & ldots \ begin{matrix} rho^{n — 2} & rho^{n — 3} \ rho^{n — 1} & rho^{n — 2} \ end{matrix} & ldots & begin{matrix} 1 & rho \ rho & 1 \ end{matrix} \ end{pmatrix}^{- 1}X} right)^{- 1} ast}X’begin{pmatrix} begin{matrix} 1 & rho \ rho & 1 \ end{matrix} & ldots & begin{matrix} rho^{n — 2} & rho^{n — 1} \ rho^{n — 3} & rho^{n — 2} \ end{matrix} \ ldots & ldots & ldots \ begin{matrix} rho^{n — 2} & rho^{n — 3} \ rho^{n — 1} & rho^{n — 2} \ end{matrix} & ldots & begin{matrix} 1 & rho \ rho & 1 \ end{matrix} \ end{pmatrix}^{- 1}y)

Обратите внимание, что дробь (frac{sigma_{u}^{2}}{1 — rho^{2}}) при расчете представленной оценки сокращается. Поэтому для вычисления оценки знать величину (sigma_{u}^{2}) не нужно.

Примечание: если коэффициент автокорреляции (rho) неизвестен, то его можно легко оценить. Например, для этого можно применить обычный МНК к исходной регрессии, получить вектор остатков и оценить регрессию ({widehat{e}}_{i} = {widehat{rho} ast e_{i — 1}}). Полученной оценки (widehat{rho}) достаточно, чтобы вычислить ОМНК-оценку вектора параметров модели. Тем самым в представленном примере для применения доступного ОМНК достаточно оценить всего один параметр ковариационной матрицы вектора оценок коэффициентов.

Пример 5.5. Гетероскедастичность и ОМНК-оценка

Рассмотрим линейную модель ({y = X}{beta + varepsilon}), для которой выполнены все предпосылки классической линейной модели множественной регрессии за одним исключением: дисперсия случайной ошибки прямо пропорциональна квадрату некоторой известной переменной

(mathit{var}{left( varepsilon_{i} right) = sigma_{i}^{2} = sigma_{0}^{2}}{z_{i}^{2} > 0}.)

(а) Найдите ковариационную матрицу вектора случайных ошибок для представленной модели.

(б) Запишите в явном виде формулу ОМНК-оценки вектора коэффициентов модели.

Решение:

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

(Omega = begin{pmatrix} begin{matrix} {sigma_{0}^{2}z_{1}^{2}} & 0 \ 0 & {sigma_{0}^{2}z_{2}^{2}} \ end{matrix} & begin{matrix} ldots & 0 \ ldots & 0 \ end{matrix} \ begin{matrix} ldots & ldots \ 0 & 0 \ end{matrix} & begin{matrix} ldots & ldots \ ldots & {sigma_{0}^{2}z_{n}^{2}} \ end{matrix} \ end{pmatrix})

(б) Обратите внимание, что при подстановке в общую формулу для ОМНК-оценки величина (sigma_{0}^{2}) сокращается, следовательно, для оценки вектора коэффициентов знать её не нужно:

({{widehat{beta}}^{} = {({X’Omega^{- 1}X})}^{- 1}}X^{‘}Omega^{- 1}{y =}left( {X’begin{pmatrix} begin{matrix} z_{1}^{2} & 0 \ 0 & z_{2}^{2} \ end{matrix} & begin{matrix} ldots & 0 \ ldots & 0 \ end{matrix} \ begin{matrix} ldots & ldots \ 0 & 0 \ end{matrix} & begin{matrix} ldots & ldots \ ldots & z_{n}^{2} \ end{matrix} \ end{pmatrix}^{- 1}X} right)^{- 1}X’begin{pmatrix} begin{matrix} z_{1}^{2} & 0 \ 0 & z_{2}^{2} \ end{matrix} & begin{matrix} ldots & 0 \ ldots & 0 \ end{matrix} \ begin{matrix} ldots & ldots \ 0 & 0 \ end{matrix} & begin{matrix} ldots & ldots \ ldots & z_{n}^{2} \ end{matrix} \ end{pmatrix}^{- 1}y)

* * *

Ещё одна важная ситуация, когда с успехом может быть применен доступный ОМНК — это модель со случайными эффектами, которую мы рассмотрим в главе, посвященной панельным данным.

Фильтр Калмана — это легко +64

Об алгоритме

Что же нам потребуется для работы фильтра Калмана?

• Нам потребуется модель системы.
• Модель должна быть линейной (об этом чуть позже).
• Нужно будет выбрать переменные, которые будут описывать состояние системы (“вектор состояния”).
• Мы будем производить измерения и вычисления в дискретные моменты времени (например, 10 раз в секунду). Нам потребуется модель наблюдения.
• Для работы фильтра достаточно данных измерений в текущий момент времени и результатов вычислений с предыдущего момента времени.

Алгоритм работает итеративно. На каждом шаге алгоритм берёт данные датчиков (с шумом и другими проблемами), вектор состояния с предыдущего шага и по этим данным оценивает состояние системы на текущем шаге. Кроме того, он еще отслеживает насколько мы можем быть уверены, что наш текущий вектор состояний соответствует истинному положению дел (разброс значений для каждой переменной в векторе).

Обычно используются следующие обозначения:

— вектор состояния;

— мера неопределенности вектора состояния. Представляет из себя ковариационную матрицу (об этом позже — это будет, наверное, самая сложная часть).

Содержимое вектора состояния зависит от фантазии разработчика и решаемой задачи. Например, мы можем отслеживать координаты объекта, а также его скорость и ускорение. В этом случае получается вектор из трёх переменных: {позиция, скорость, ускорение} (для одномерного случая; для 3D мира будет по одному такому набору для каждой оси, то есть, 9 значений в векторе)

Математическая модель системы / процесса

Мы имеем дело с динамической системой, т.е. состояние системы меняется со временем. Имея модель системы, фильтр Калмана может предугадывать, каким будет состояние системы в следующий момент времени. Именно это позволяет фильтру так эффективно устранять шум и оценивать параметры, которые не наблюдаются (не измеряются) напрямую.

Фильтр Калмана накладывает ограничения на используемые модели: это должны быть дискретные модели в пространстве состояний. А ещё они должны быть линейными.

Такую модель удобно представлять в виде разностного матричного уравнения:

Давайте разберём это уравнение подробно. В первую очередь, нас интересует первое слагаемое (

) — это как раз модель эволюции процесса. А матрица

(также встречаются обозначения

,

) — называется матрицей процесса (state transition matrix). Она задаёт систему линейных уравнений, описывающих, как получается новое состояние системы из предыдущего.

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

Первая строка матрицы — хорошо знакомое уравнение

. Аналогично, вторая строка матрицы описывает изменение скорости. Третья строка описывает изменение ускорения.

Модель наблюдения

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

— это вектор измерения/наблюдения.Это значения, получаемые с датчиков системы.

Первое слагаемое

— модель, связывающая вектор состояния

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

из

, а эта модель получает

из

. Но это действительно так. В частности, это необходимо потому, что некоторые переменные состояния из

могут отсутствовать в

).

Второе слагаемое

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

— ковариационная матрица, соответствующая вектору

.

Вернёмся к нашему примеру. Пусть у нас на роботе установлен один единственный датчик — GPS приёмник (“измеряет” положение). В этом случае матрица

будет выглядеть следующим образом:

Строки матрицы соответствуют переменным в векторе состояния, столбцы — элементам вектора измерений. В первой строке матрицы находится значение “1” так как единица измерения положения в векторе состояния совпадает с единицей измерения значения в векторе измерений. Остальные строки содержат “0” потому что переменные состояния соответствующие этим строкам не измеряются датчиком.

Ковариационные матрицы и где они обитают

Для настройки фильтра нам потребуется заполнить несколько ковариационных матриц:

,

и

.

На каждом шаге фильтр Калмана строит предположение о состоянии системы, исходя из предыдущей оценки состояния и данных измерений. Если неопределенности вектора состояния выше, чем ошибка измерения, то фильтр будет выбирать значения ближе к данным измерений. Если ошибка измерения больше оценки неопределенности состояния, то фильтр будет больше “доверять” данным моделирования. Именно поэтому важно правильно подобрать значения ковариационных матриц — основного инструмента настройки фильтра.

Рассмотрим каждую матрицу подробнее:

— ковариационная матрица состояния

Квадратная матрица, порядок матрицы равен размеру вектора состояния

Как уже было сказано выше, эта матрица определяет “уверенность” фильтра в оценке переменных состояния. Алгоритм самостоятельно обновляет эту матрицу в процессе работы. Однако нам нужно установить начальное состояние, вместе с исходным предположением о векторе состояния.

Во многих случаях нам неизвестны значения ковариации между переменными для изначального состояния (элементы матрицы, расположенные вне главной диагонали). Поэтому можно проигнорировать их, установив равными 0. Фильтр самостоятельно обновит значения в процессе работы. Если же значения ковариации известны, то, конечно же, стоит использовать их.

Дисперсию же проигнорировать не выйдет. Необходимо установить значения дисперсии в зависимости от нашей уверенности в исходном векторе состояния. Для этого можно воспользоваться правилом трёх сигм: значение случайной величины попадает в диапазон

с вероятностью 99.7%.

— ковариационная матрица шума измерений

Квадратная матрица, порядок матрицы равен размеру вектора наблюдения (количеству измеряемых параметров).

Во многих случаях можно считать, что измерения не коррелируют друг с другом. В этом случае матрица

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

— ковариационная матрица ошибки модели

Квадратная матрица, порядок матрицы равен размеру вектора состояния.

С этой матрицей обычно возникает наибольшее количество вопросов. Что означает ошибка модели? Каков смысл этой матрицы и за что она отвечает? Как заполнять эту матрицу? Рассмотрим всё по порядку.

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

Если установить очень маленькое значение

, то этап предсказания будет слабо увеличивать неопределенность оценки. Это означает, что мы считаем, что наша модель точно описывает процесс.

Если же установить большое значение

, то этап предсказания будет сильно увеличивать неопределенность оценки. Таким образом, мы показываем что модель может содержать неточности или неучтенные факторы.

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

Допустим, мы отслеживаем перемещение робота, используя модель равноускоренного движения, и вектор состояния содержит следующие переменные: положение x, скорость v и ускорение a. Однако, наша модель не учитывает, что на дороге встречаются неровности.

Когда робот проходит неровность, показания датчиков и предсказание модели начнут расходиться. Структура матрицы

будет определять, как фильтр отреагирует на это расхождение.

Мы можем выдвинуть различные предположения относительно природы шума. Для нашего примера с равноускоренным движением логично было бы предположить, что неучтённые факторы (неровность дороги) в первую очередь влияют на ускорение. Этот подход применим ко многим структурам модели, где в векторе состояния присутствует переменная и несколько её производных по времени (например, положение и производные: скорость и ускорение). Матрица

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

Простейший подход

В некоторых случаях прибегают к грубому упрощению: устанавливают все элементы матрицы

равными нулю, за исключением элементов, соответствующих максимальным порядкам производных переменных состояния.

Действительно, если рассчитать

по одному из приведённых выше методов, при достаточно малых значениях

, значения элементов матрицы

оказываются очень близкими к нулю.

Т.е. для нашей модели равноускоренного движения можно взять матрицу

следующего вида:

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

таким образом для любых важных задач без весомых причин.

Важное замечание

Во всех примерах выше используется вектор состояния

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

Рассмотрим вектор состояния

Матрица

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

и

. Остальные элементы матрицы будут равны нулю.

Дисперсия, соответствующая наивысшим порядкам производных

и

, будет находиться на 3-ей и 5-ой позициях на главной диагонали матрицы.

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

Пример кода

Нет смысла изобретать велосипед и писать свою собственную реализацию фильтра Калмана, когда существует множество готовых библиотек. Я выбрал язык python и библиотеку filterpy для примера.

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

Моделирование данных датчиков

Начнём с равномерного движения:

import numpy as np
import numpy.random

# Моделирование данных датчика
def simulateSensor(samplesCount, noiseSigma, dt):
   # Шум с нормальным распределением. мат. ожидание = 0, среднеквадратичное отклонение = noiseSigma
   noise = numpy.random.normal(loc = 0.0, scale = noiseSigma, size = samplesCount)

   trajectory = np.zeros((3, samplesCount))

   position = 0
   velocity = 1.0
   acceleration = 0.0

   for i in range(1, samplesCount):
       position = position + velocity * dt + (acceleration * dt ** 2) / 2.0
       velocity = velocity + acceleration * dt
       acceleration = acceleration

       trajectory[0][i] = position
       trajectory[1][i] = velocity
       trajectory[2][i] = acceleration

   measurement = trajectory[0] + noise

   return trajectory, measurement # Истинное значение и данные "датчика" с шумом

Визуализируем результаты моделирования:

import matplotlib.pyplot as plt
dt = 0.01
measurementSigma = 0.5
trajectory, measurement = simulateSensor(1000, measurementSigma, dt)

plt.title("Данные датчика")
plt.plot(measurement, label="Измерение", color="#99AAFF")
plt.plot(trajectory[0], label="Истинное значение", color="#FF6633")
plt.legend()
plt.show()

Реализация фильтра

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

У нас единственный датчик, который напрямую измеряет положение. Поэтому модель наблюдения получается очень простой:

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

Однако, мы не уверены, что это именно так. Поэтому установим матрицу ковариации для начального состояния с большими значениями на главной диагонали:

Я воспользовался функцией библиотеки filterpy для расчёта ковариационной матрицы ошибки модели: filterpy.common.Q_discrete_white_noise. Эта функция использует модель непрерывного белого шума.

import filterpy.kalman
import filterpy.common
import matplotlib.pyplot as plt
import numpy as np
import numpy.random
from Simulator import simulateSensor # моделирование датчиков

dt = 0.01                       # Шаг времени
measurementSigma = 0.5          # Среднеквадратичное отклонение датчика
processNoise = 1e-4             # Погрешность модели

# Моделирование данных датчиков
trajectory, measurement = simulateSensor(1000, measurementSigma, dt)

# Создаём объект KalmanFilter
filter = filterpy.kalman.KalmanFilter(dim_x=3,      # Размер вектора стостояния
                                     dim_z=1)      # Размер вектора измерений

# F - матрица процесса - размер dim_x на dim_x - 3х3
filter.F = np.array([ [1,   dt,     (dt**2)/2],
                     [0,   1.0,    dt],
                     [0,   0,      1.0]])


# Матрица наблюдения - dim_z на dim_x - 1x3
filter.H = np.array([[1.0, 0.0, 0.0]])

# Ковариационная матрица ошибки модели
filter.Q = filterpy.common.Q_discrete_white_noise(dim=3, dt=dt, var=processNoiseVariance)

# Ковариационная матрица ошибки измерения - 1х1
filter.R = np.array([[measurementSigma*measurementSigma]])

# Начальное состояние.
filter.x = np.array([0.0, 0.0, 0.0])

# Ковариационная матрица для начального состояния
filter.P = np.array([[10.0, 0.0,  0.0],
                    [0.0,  10.0, 0.0],
                    [0.0,  0.0,  10.0]])

filteredState = []
stateCovarianceHistory = []

# Обработка данных
for i in range(0, len(measurement)):
   z = [ measurement[i] ]                      # Вектор измерений
   filter.predict()                            # Этап предсказания
   filter.update(z)                            # Этап коррекции

   filteredState.append(filter.x)
   stateCovarianceHistory.append(filter.P)

filteredState = np.array(filteredState)
stateCovarianceHistory = np.array(stateCovarianceHistory)

# Визуализация
plt.title("Kalman filter (3rd order)")
plt.plot(measurement, label="Измерение", color="#99AAFF")
plt.plot(trajectory[0], label="Истинное значение", color="#FF6633")
plt.plot(filteredState[:, 0], label="Оценка фильтра", color="#224411")
plt.legend()
plt.show()

Бонус — сравнение различных порядков моделей

Сравним поведение фильтра с моделями разного порядка. Для начала, смоделируем более сложный сценарий поведения робота. Пусть робот находится в покое первые 20% времени, затем движется равномерно, а затем начинает двигаться равноускоренно:

# Моделирование данных датчика
def simulateSensor(samplesCount, noiseSigma, dt):
   # Шум с нормальным распределением. мат. ожидание = 0, среднеквадратичное отклонение = noiseSigma
   noise = numpy.random.normal(loc = 0.0, scale = noiseSigma, size = samplesCount)

   trajectory = np.zeros((3, samplesCount))

   position = 0
   velocity = 0.0
   acceleration = 0.0

   for i in range(1, samplesCount):
       position = position + velocity * dt + (acceleration * dt ** 2) / 2.0
       velocity = velocity + acceleration * dt
       acceleration = acceleration

       # Переход на равномерное движение
       if(i == (int)(samplesCount * 0.2)):
           velocity = 10.0

       # Переход на равноускоренное движение
       if (i == (int)(samplesCount * 0.6)):
           acceleration = 10.0

       trajectory[0][i] = position
       trajectory[1][i] = velocity
       trajectory[2][i] = acceleration

   measurement = trajectory[0] + noise

   return trajectory, measurement # Истинное значение и данные "датчика" с шумом

В предыдущем примере мы использовали модель, содержащую переменную (положение) и две производных её по времени (скорость и ускорение). Посмотрим, что будет, если избавиться от одной или обеих производных:

# Создаём объект KalmanFilter
filter = filterpy.kalman.KalmanFilter(dim_x=2,      # Размер вектора стостояния
                                     dim_z=1)      # Размер вектора измерений

# F - матрица процесса - размер dim_x на dim_x - 2х2
filter.F = np.array([ [1,   dt],
                     [0,   1.0]])


# Матрица наблюдения - dim_z на dim_x - 1x2
filter.H = np.array([[1.0, 0.0]])

filter.Q = [[dt**2,       dt],
           [   dt,      1.0]] * processNoise

# Начальное состояние.
filter.x = np.array([0.0, 0.0])

# Ковариационная матрица для начального состояния
filter.P = np.array([[8.0, 0.0],
                    [0.0, 8.0]])

# Создаём объект KalmanFilter
filter = filterpy.kalman.KalmanFilter(dim_x=1,      # Размер вектора стостояния
                                     dim_z=1)      # Размер вектора измерений

# F - матрица процесса - размер dim_x на dim_x - 1х1
filter.F = np.array([ [1.0]])


# Матрица наблюдения - dim_z на dim_x - 1x1
filter.H = np.array([[1.0]])

# Ковариационная матрица ошибки модели
filter.Q = processNoise 

# Ковариационная матрица ошибки измерения - 1х1
filter.R = np.array([[measurementSigma*measurementSigma]])

# Начальное состояние.
filter.x = np.array([0.0])

# Ковариационная матрица для начального состояния
filter.P = np.array([[8.0]])

Сравним результаты:

На графиках сразу заметно, что модель первого порядка начинает отставать от истинного значения на участках равномерного движения и равноускоренного движения. Модель второго порядка успешно справляется с участком равномерного движения, но так же начинает отставать на участке равноускоренного движения. Модель третьего порядка справляется со всеми тремя участками.

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

Нелинейные модели и фильтр Калмана

Почему фильтр Калмана не работает для нелинейных моделей и что делать

Всё дело в нормальном распределении. При применении линейных преобразованийк нормально распределенной случайной величине, результирующее распределение будет представлять собой нормальное распределение, или будет пропорциональным нормальному распределению. Именно на этом принципе и строится математика фильтра Калмана.

Есть несколько модификаций алгоритма, которые позволяют работать с нелинейными моделями.

Например:

Extended Kalman Filter (EKF) — расширенный фильтр Калмана. Этот подход строит линейное приближение модели на каждом шаге. Для этого требуется рассчитать матрицу вторых частных производных функции модели, что бывает весьма непросто. В некоторых случаях, аналитическое решение найти сложно или невозможно, и поэтому используют численные методы.

Unscented Kalman Filter (UKF). Этот подход строит приближение распределения получающегося после нелинейного преобразования при помощи сигма-точек. Преимуществом этого метода является то, что он не требует вычисления производных.

Мы рассмотрим именно Unscented Kalman Filter

Unscented Kalman Filter и почему он без запаха

Основная магия этого алгоритма заключается в методе, который строит приближение распределения плотности вероятности случайной величины после прохождения через нелинейное преобразование. Этот метод называется unscented transform — сложнопереводимое на русский язык название. Автор этого метода, Джеффри Ульман, не хотел, чтобы его разработку называли “Фильтр Ульмана”. Согласно интервью, он решил назвать так свой метод после того как увидел дезодорант без запаха (“unscented deodorant”) на столе в лаборатории, где он работал.

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

Для использования UKF не придётся реализовывать какие-либо дополнительные вычисления, за исключением моделей системы. В общем виде, нелинейная модель не может быть представлена в виде матрицы, поэтому мы заменяем матрицы

и

на функции

и

. Однако смысл этих моделей остаётся тем же.

Реализуем unscented Kalman filter для линейной модели из прошлого примера:

import filterpy.kalman
import filterpy.common
import matplotlib.pyplot as plt
import numpy as np
import numpy.random
from Simulator import simulateSensor, CovarianceQ

dt = 0.01
measurementSigma = 0.5
processNoiseVariance = 1e-4

# Функция наблюдения - аналог матрицы наблюдения
# Преобразует вектор состояния x в вектор измерений z
def measurementFunction(x):
   return np.array([x[0]])

# Функция процесса - аналог матрицы процесса
def stateTransitionFunction(x, dt):
   newState = np.zeros(3)
   newState[0] = x[0] + dt * x[1] + ( (dt**2)/2 ) * x[2]
   newState[1] = x[1] + dt * x[2]
   newState[2] = x[2]

   return newState


trajectory, measurement = simulateSensor(1000, measurementSigma)

# Для unscented kalman filter необходимо выбрать алгоритм выбора сигма-точек
points = filterpy.kalman.JulierSigmaPoints(3, kappa=0)

# Создаём объект UnscentedKalmanFilter
filter = filterpy.kalman.UnscentedKalmanFilter(dim_x = 3,
                                              dim_z = 1,
                                              dt = dt,
                                              hx = measurementFunction,
                                              fx = stateTransitionFunction,
                                              points = points)



# Ковариационная матрица ошибки модели
filter.Q = filterpy.common.Q_discrete_white_noise(dim=3, dt=dt, var=processNoiseVariance)

# Ковариационная матрица ошибки измерения - 1х1
filter.R = np.array([[measurementSigma*measurementSigma]])

# Начальное состояние.
filter.x = np.array([0.0, 0.0, 0.0])

# Ковариационная матрица для начального состояния
filter.P = np.array([[10.0, 0.0,  0.0],
                    [0.0,  10.0, 0.0],
                    [0.0,  0.0,  10.0]])

filteredState = []
stateCovarianceHistory = []

for i in range(0, len(measurement)):
   z = [ measurement[i] ]
   filter.predict()
   filter.update(z)

   filteredState.append(filter.x)
   stateCovarianceHistory.append(filter.P)

filteredState = np.array(filteredState)
stateCovarianceHistory = np.array(stateCovarianceHistory)

plt.title("Unscented Kalman filter")
plt.plot(measurement, label="Измерение", color="#99AAFF")
plt.plot(trajectory[0], label="Истинное значение", color="#FF6633")
plt.plot(filteredState[:, 0], label="Оценка фильтра", color="#224411")
plt.legend()
plt.show()

Разница в коде минимальна. Мы заменили матрицы F и H на функции f(x) и h(x). Это позволяет использовать нелинейные модели системы и/или наблюдения:

# Функция наблюдения - аналог матрицы наблюдения
# Преобразует вектор состояния x в вектор измерений z
def measurementFunction(x):
   return np.array([x[0]])

# Функция процесса - аналог матрицы процесса
def stateTransitionFunction(x, dt):
   newState = np.zeros(3)
   newState[0] = x[0] + dt * x[1] + ( (dt**2)/2 ) * x[2]
   newState[1] = x[1] + dt * x[2]
   newState[2] = x[2]

   return newState

Также, появилась строчка, устанавливающая алгоритм генерации сигма-точек

points = filterpy.kalman.JulierSigmaPoints(3, kappa=0)

Этот алгоритм определяет точность оценки распределения вероятности при прохождении через нелинейное преобразование. К сожалению, существуют только общие рекомендации относительно генерации сигма-точек. Поэтому для каждой отдельной задачи значения параметров алгоритма подбираются экспериментальным путём.

Ожидаемый результат — график оценки положения практически не отличается от обычного фильтра Калмана.

В этом примере используется линейная модель. Однако мы могли бы использовать нелинейные функции. Например, мы могли бы использовать следующую реализацию:

g = 9.8
# Вектор состояния - угол наклона
# Вектор измерений - ускорение вдоль осей X и Y
def measurementFunction(x):
  measurement = np.zeros(2)
  measurement[0] = math.sin(x[0]) * g
  measurement[1] = math.cos(x[0]) * g
  return measurement

Такую модель измерений было бы невозможно использовать в случае с линейным фильтром Калмана

Вместо заключения

За рамками статьи остались теоретические основы фильтра Калмана. Однако объем материала по этой теме ошеломляет. Сложно выбрать хороший источник. Я бы хотел рекомендовать замечательную книгу от автора библиотеки filterpy Roger Labbe (на английском языке). В ней доступно описаны принципы работы алгоритма и сопутствующая теория. Книга представляет собой документ Jupyter notebook, который позволяет в интерактивном режиме экспериментировать с примерами.

Литература

> Roger Labbe — Kalman and Bayesian Filters in Python
> Wikipedia

В течение последних лет методы
нестационарной (во времени) фильтрации и сглаживания широко использовались в
различных задачах. Задачи фильтрации и сглаживания подразумевают оценивание сигнала
по наблюдениям, состоящим из аддитивной смеси сигнала и шума. Различие этих
двух задач состоит в моментах времени оценки сигнала и в интервалах времени
наблюдения. В задачах фильтрации оценку получают в конце интервала наблюдения,
при сглаживании сигнал оценивается в некоторые моменты внутри интервала
наблюдения. В гл. 7 и § 8.2, 8 3 эти проблемы основательно обсуждены.

Две трудности
препятствуют практическому использованию уравнений фильтрации и сглаживания –
это выбор математической модели вектора состояний оцениваемого сигнала и выбор
соответствующих ковариационных матриц сигнала и наблюдения, используемых при
решении этих задач. Желательно, чтобы модель сигнала описывала его достаточно
полно и в то же время была проста настолько, чтобы алгоритмы фильтрации
поддавались численному решению Матрицы ковариаций для сигнала и шума измерений
также желательно иметь такие, чтобы они соответствовали реальным процессам. К сожалению,
полное статистическое описание подобных процессов, как правило, невозможно.

Основным показателем расхождения
(модели и сигнала и матриц ковариаций является анализ чувствительности
ковариационной матрицы ошибок оценивания. Анализ чувствительности в задачах
фильтрации при больших ошибках был рассмотрен для случая дискретных наблюдений
в [60], [87], [173], [174], а для непрерывного времени в [175] Детальный анализ
чувствительности в задачах сглаживания дан в [82], [84] В этом разделе
построим алгоритмы, необходимые для анализа чувствительности в задачах
сглаживания и фильтрации. Эти алгоритмы будут получены для задач сглаживания в
фиксированной точке и сглаживания на временном фиксированном интервале. Отсюда
как частные случаи будут получены алгоритмы чувствительности задач фильтрации.

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

Предположим,
что процесс описывается моделью

                                                                   (8.178)

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

                                                                              (8.179)

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

                                                        (8.180)

может быть описан уравнениями,
содержащимися в табл 7.4.

                                (8.181)

                                                       (8.182)

  (8.182 а)

где  — оптимальная оценка вектора состояний в
момент t при предполагаемых параметрах .

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

                                                               (8.183)

дано в табл. 8.3. Оценка вектора
состояний удовлетворяет уравнению

.           (8.184)

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

                              (8.185)

Начальные условия
для (8.185)
вытекают из решения задачи фильтрации. Как уже известно, ур-ния (8.184) и (8 185) решаются в обратном времени
для всех  и
решением является оценка сглаживания, а также соответствующая ковариационная
матрица ошибок В ажио отметить, что для того чтобы найти оценку сглаживания,
нет необходимости вычислять (ковариационную’ матрицу ошибок.

Указанные
алгоритмы оптимальной фильтрации и сглаживания обеспечивают минимальный риск
только в том случае, если модель совпадает с реальным процессом. Если же этого
не происходит, то матрица
ковариаций ошибок оценивания уже не является истинной (ковариационной матрицей
ошибок Таким образом, действительная ковариационная матрица ошибок может
выявлять ошибки, возникающие в результате использования некорректной модели.
Поэтому эта матрица играет первостепенную роль в любом анализе ошибок и
исследовании чувствительности алгоритмов фильтрации и сглаживания Как сейчас
будет показано, для действительной матрицы ковариаций можно получить
специальное дифференциальное уравнение

Сначала предположим,
что реальный физический процесс точно описывается уравнением

                                          (8.187)

а вектор измерений равен

                                                (8.188)

и величины  и  являются точными. Теперь
найдем уравнения для действительной ошибки фильтрации

                                                          (8.189)

и ошибки сглаживания

                                                     (8.190)

Используя
введенные соотношения, можно получить дифференциальные уравнения, которым
должны удовлетворять ошибки оценивания. В случае фильтрации имеем

      (8.191)

а для сглаживания

          (8.192)

где

                                 (8.193)

В ур-ниях (8 191) и (8.192) для
упрощения записи у некоторых функций не указан аргумент (время) Образовывая
новый вектор состояний, можно векторные дифференциальные уравнения для
действительной ошибки фильтрации сглаживания и ошибки модели свести в одно
дифференциальное уравнении. Этот вектор можно определить как

,                                                          (8.194)

а новый входной вектор

                .                                                             (8.195)

Тогда новое уравнение для вектора
состояний можно записать как

                             (8.196)

        (8.197)

                                                        (8.198)

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

Дифференциальное
уравнение для  имеет
следующее решение

                 (8.199)

где  — переходная матрица состояний (8.199). Используя
(8.196) и (8.199), а также свойство дельта-функции, имеем

          (8.200)

                                                                     (8.201)

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

Расширяя матрицу ковариаций для получим

                                          (8.202)

где  В алгоритмах фильтрации действительная
ковариационная матрица ошибок оценивания определяется решением
дифференциального матричного уравнения

                          (8.203)

где  удовлетворяет условию

                        (8.204)

а  условию

                                              (8.205) 

Эти соотношения легко получить
методами, развитыми в § 3.5. Начальными условиями для предыдущих уравнений
являются

                                                                              (8.206)

Если матрицы  и  одинаковы для модели и
действительного физического процесса, то ур-ние (8.203) приобретает такой вид,
что для его решения не нужны ур-ния (8.204) и (8.205). Наконец, если  и  то (8.203) переходит в
(8.183), являющееся строгим уравнением для ковариационной матрицы ошибок, если
модель и априорные данные верны.

Дополнительные дифференциальные
уравнения:

                               (8.207)

              (8.208)

                               (8.209)

                            (8.210)

                     (8.211)

с начальными условиями:

                                               (8.212)

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

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

                                                   (8.213)

                                                              (8.214)

Поскольку правые части этих тождеств удовлетворяют
дифференциальным уравнениям и начальным условиям левой часть, то указанные
тождества справедливы. Использование этих тождеств в ур-нии (8.207) дает

              (8.215)

что является эквивалентным выражением для полученных
ранее соотношений (8.207) — (8.209). Для нахождения ковариационной матрицы
сглаживания ур-ния (8 210), (8.211) и (8.215) должны решаться в обратном
времени для всех .
Если модель действительно отражает реальный физический процесс, то ур-ние
(8.215) переходит в (8.186).

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

                 (8.216)

для случаев фильтрации и сглаживания соответственно;  — разность между значением параметра модели,
используемой для оценивания, и параметра действительного процесса. Величины  и  вычисляются при действительном значении параметра , a  и  — при
предполагаемом его значении .

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

                                                (8.217)

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

,                                                                 (8.218)

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

Если взять
частную производную от обеих частей ур-ния (8.203), поменять порядок
дифференцирования и множество параметров действительного процесса заменить на
множество параметров модели фильтрации, то получим дифференциальное уравнение
для чувствительности в малом, которое имеет вид

                       (8.219)

где

   (8.220)           

Начальным условием для этого выражения является
соотношение

                                                                (8.221) 

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

      Матричная функция
чувствительности
в малом позволяет найти еще две величины,
представляющие интерес. Можно найти матрицу разности между действительной
ковариационной матрицей и ковариационной матрицей моде

                                                (8.222) 

где  — число компонент вектора , a  —
дифференциальный оператор. Матрицу чувствительности в малом можно использовать
также для нахождения чувствительности в малом функции риска тем же методом,
который использовался при анализе чувствительности в большом. В результате
чувствительность функции риска равна

                                                                       (8.223) 

Полученное выше выражение для
чувствительности функции риска имеет такой же вид, который использовался в
[202].

Для задач сглаживания
также можно определить матричную функцию чувствительности

                                                                             (8.224) 

Так же, как и в случае фильтрации,
«можно получить для этой функции следующие дифференциальные уравнения:

          (8.225) 

                             (8.226)

где

                        (8.227) 

Кроме того,

                                          (8.228) 

Начальными условиями являются:

                                                                              (8.229) 

Аналогичные результаты для анализа
чувствительности в малом можно получить и для алгоритмов сглаживания на
фикшрован ном интервале. В табл. 8.9 ,и 8.10 приведены алгоритмы анализа ошибок
и чувствительности для задач фильтрации и сглаживания на фиксированном
интервале в непрерывном времени. Используя ту же методику, можно получить и
дискретный аналог выведенных алгоритмов. Без вывода они даны в табл. 8.11 и
8.12.

Таблица 8.9. Алгоритм
анализа ошибок и чувствительности для фильтрации в непрерывном времени

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

Таблица 8.11. Дискретные
алгоритмы анализа ошибок и чувствительности фильтрации

Таблица 8.12.  Дискретные алгоритмы анализа ошибок
и чувствительности для сглаживания на фиксированном интервале (алгоритмы
фильтрации табл. 8 11 составляют общую часть алгоритмов этой таблицы)

Поскольку получение предыдущих результатов заняло много места и времени,
мы рассмотрим один простой пример, чтобы продемонстрировать разработанные идеи.
Более сложные примеры можно найти в [82, 84].

Пример 8.3. Рассмотрим
простой пример, в котором предположим, что модель процесса
описывается соотношениями:

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

Уравнение для ковариационной матрицы ошибок фильтрации имеет вид  Решение этого
уравнения  и,
следовательно, коэффициент усиления .

Если
предположить, что параметры  и  являются
действительными, то ур-ние (8.203) для действительной ковариационной матрицы
ошибок примет вид

где  удовлетворяет уравнению

и  . Ho  есть
дисперсия состояния и . Начальные условия . Решение приведенных
уравнений дает:

Непосредственное
использование (8.216) дает выражение для чувствительности в большом алгоритмов
фильтрации:

Дифференциальные уравнения
матричных функций чувствительности в малом можно получить из (8.219) н (8.220)
в виде

где  соответствует ,  или ,
при этом

с начальными условиями . Эти дифференциальные
уравнения первого порядка легко решаются:

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

На рис. 8.1 и 8.2
приведены для данного примера функции чувствительности в большом и в малом.
Здесь выражения для чувствительностей в большом и малом одинаковы для исходной
и вычисленной ковариационной матриц.

Рис. 8.1. Функции
чувствительности дисперсии ошибки в зависимости от изменений начальной дисперсии
ошибки и дисперсии измерений: фильтрация;  ______ — сглаживание;

Рис. 8.2. Функции
чувствительности при изменении константы измерений:

Зависимость чувствительности от параметра  приведена на рис. 8.3,
иллюстрирующем различие между чувствительностью в малом и большом.

В этом примере ковариационная
матрица ошибок сглаживания  при конечном условии . Отсюда следует, что
ковариационная матрица ошибок не зависит от времени и получается из (8.203) для
момента времени .
Значения чувствительностей также постоянны и равны чувствительностям для
фильтра в момент .
Пунктирные кривые на рис. 8,1 представляют две функции чувствительности для
задачи сглаживания.

Рис.
8.3. Сравнение чувствительности в большом и малом при изменении константы
измерений.

«Вычисляемая»
дисперсия ошибки фильтрации всегда уменьшается при увеличении времени
наблюдения. Для численного анализа выражений для действительной дисперсии
положим .
Тогда «вычисленная» дисперсия ошибки равна , а следовательно, «вычисленная» дисперсия
сглаживания  и,
как и следовало ожидать, от  не зависит. В этом случае действительная
дисперсия ошибки фильтрации определится выражением . Если , то это выражение совпадает с
выражением для минимума дисперсии ошибки .

     Интересно
отметить, что действительная ошибка (дисперсия) не достигает нуля при больших t до тех пор, пока  не станет равным
нулю, что означает отсутствие ошибки в модели наблюдения. Для больших  приближенно можно
считать, что действительная дисперсия ошибки равна . На рис. 8.4 приведены графики
действительных дисперсий ошибок и «вычисленных» дисперсий для случая, когда
априорная статистика задана верно .

       Если
модель наблюдения не содержит ошибок, , то дисперсия ошибки . Если отношение , то имеем минимум дисперсии
ошибки, соответствующий .

Анализ
ошибок и чувствительности линейных алгоритмов сглаживания в фиксированной
точке.

Предположим,
что непрерывный процесс может быть задан векторным дифференциальным     уравнением
(8.178), а модель наблюдения линейна [см. (8.179)]

Рис. 8.4. Дисперсия ошибки фильтрации при разных ,;

Оценка
сглаживания в фиксированной точке удовлетворяет дифференциальному уравнению

                          (8.230)
   

где

,                                           (8.231)                                                                                                                  

а   — матрица (обозначенная
так для удобства), которая удовлетворяет уравнению

                         (8.232)

Начальным
условием служит оценка  фильтрации в момент  и

.                                                              (8.233)

Ковариационная
матрица ошибки оценивания может быть получена из уравнения

                            (8.234)

Эта матрица ошибок удовлетворяет
ур-нию (8.234) только в том случае,
если модель и ее априорные характеристики совпадают с действительным процессом
и его характеристиками. На практике обычно это условие не выполняется.
Поскольку действительная ковариационная матрица ошибок оценивания позволяет
выявить ошибки при использовании неадэкватной модели, она играет основную роль в анализе ошибок и чувствительности алгоритмов сглаживания

В непрерывном времени
действительная ковариационная матрица вычисляется следующим образом. Сначала
постулируем модель для действительного процесса в том же виде, что и в ур-ниях
(8.178), (8.179), за исключением того, что не ставим черты сверху. Далее
находим дифференциальное уравнение для действительной ошибки оценивания
подстановкой соотношений:

                     (8.235)

в (8.230), где
действительные ошибки оценок фильтрации и сглаживания определены выражениями:

                 (8.236)

Тогда результирующее
дифференциальное уравнение для ошибки оценивания имеет вид

                               (8.237)

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

,                                                             (8.238)

где

                                          (8.239)

Матрицы   и  были определены в предыдущем
параграфе. Далее выпишем решение ур-ния (8.238):

,                    (8.240)

где  удовлетворяет матричному
дифференциальному уравнению

.                                        (8.241)

Затем необходимо усреднить
произведение вектора  на транспонированный вектор . Тогда получим

  (8.242)

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

                                       (8.243)

Если взять внутренний
интеграл в последнем слагаемом (8.242), то, используя свойства дельта-функции,
можно получить

        (8.244)

Члены этого уравнения
содержат искомые ковариационные матрицы. Однако численное решение этого
уравнения затруднительно. Эти трудности можно обойти, если найти дифференциальное
уравнение для .
Его можно получить следующим образом. Сначала берется производная от обеих
частей ур-ния (8.244). Далее в полученное уравнение подставляется ур-ние
(8.241). Наконец используя (8.240), находим

                              (8.245)

По определению
имеем

                                (8.246)

где

                  (8.247)

а  определено соотношением (8.204).
Дифференциальное уравнение для действительной ковариационной матрицы ошибки
сглаживания имеет вид

             (8.248)

где  и  удовлетворяют уравнениям:

    (8.249)

                                           (8.250)

Начальными условиями являются:

                                                     (8.251)

Дифференциальные уравнения для ,  и  даны в предыдущем разделе
совместно с результатами фильтрации.

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

                                                           (8.252)

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

Дифференциальная матричная
функция чувствительности (чувствительность в малом) ‘находится из выражения

                                                                       (8.253)

Символом  обозначено множество
параметров действительного процесса;  — соответствующее множество параметров
модели оценивания;  -й элемент множества . Если разности между
истинными параметрами и параметрами модели малы, то функции  чувствительности в
малом могут быть использованы для нахождения матрицы разностей между
ковариационными матрицами модели и процесса. Эта матрица определяется как  где  — число
параметров множества , а — дифференциал -го параметра множества .

      Прямое вычисление (8.254) если
не невозможно, то затруднительно, потому что почти не существует аналитических
решений для действительных матриц ковариаций. Если взять частную производную от
обеих частей ур-ния (8.248), то можно найти дифференциальные уравнения,
решением которых является искомая функция чувствительности. Если после этого поменять
порядок дифференцирования, то получим

   (8.254)

где  и  удовлетворяют уравнениям:

              (8.255)

                                                 (8.256)                                                                                            

В ур-нии (8 255) член  появляется в
результате анализа чувствительности в малом задачи фильтрации, а ур-ние (8.219)
представляет алгоритм для его вычисления. Для указанных трех уравнений
начальными условиями являются:

                   (2.857)

Пример 8.4. Рассмотрим упрощенную
модель, в которой принимаемый амплитудно модулированный сигнал аддитивно
взаимодействует с белым шумом. Пусть полезный сигнал имеет вид где  и  — нормальные случайные величины
с нулевыми средними и дисперсиями  и  соответственно. Этот сигнал модулирует
несущую, так что наблюдаемое колебание имеет вид  Здесь предполагается, что  — белый нормальный шум с нулевым
средним и интенсивностью . Необходимо оценить сигнал в момент t по наблюдаемой реализации на
интервале от нуля до . Задачу можно решить методом сглаживания
в фиксированной точке, формулируя ее в терминах состояний. Такое
переопределение задачи приводит в конечном счете к векторному дифференциальному
уравнению.

где  элемент  есть сигнал . Вектор измерений
равен   

На рис. 8.5-8.7 представлены
результаты некоторых вычислений. Здесь дисперсии величин  и  равны
единице, а дисперсия шума равна 0,25. Частоты сигнала и несущей равны  6,28 и
62,8 1/с соответственно. В этих предположениях дисперсия ошибок иллюстрируется
рис. 8.5.

Рис. 8.5. Дисперсия
ошибок оценки сигнала с момента 0,2 с для различных значений дисперсии шума
наблюдений: 1 — действительная дисперсия; 2 — предполагаемая дисперсия

Рис. 8.6.
Чувствительность дисперсии ошибки оценки сигнала к 10%-ному изменению несущей
частоты в зависимости от времени наблюдения

Рис. 8 7. Дисперсия
ошибки оценки сигнала для различных значений несущей частоты: 1— действительная
дисперсия; 2 — предполагаемая дисперсия

Здесь дисперсия ошибки предсказания
сигнала в момент 0,2 с [вычисленная по (8.248)] приведена как функция времени,
при котором появляются иовые наблюдения. На рис 8 6 представлен график функции
чувствительности в большом дисперсии ошибки оценки для случая, когда истинная
частота сигнала на 10% выше предполагаемой. Влияние ошибки в выборе частоты
несущей показано на рис. 8.7.

          Алгоритмы анализа для
сглаживания в фиксированной точке в непрерывном времени сведены в табл. 8ЛЗ, а
для дискретного случая в табл. 8.14.

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