Идентификация объектов и систем управления .Теория оценивания

Идентификация объектов и систем управления Теория оценивания Андрей Масленников Москва, 2018 Тема лекций: Теория оценивания 1. Теория оценивания. 2. Непараметрическое оценивание. 3. Параметрическое оценивание. 3.1. Точечные оценки. 3.2. Интервальные оценки. 3.3. Метод минимизации среднеквадратичной ошибки. 3.4. Метод наименьших квадратов. 3.5. Метод максимального подобия. 3.6. Метод апостериорного максимума. 4. Информационные критерии и выбор модели. 5. Фильтры Калмана. 5.1. Линейный фильтр Калмана. 5.2. Расширенный фильтр Калмана. 5.3. Сигма-точечный фильтр Калмана. 1 / 90 Теория оценивания Теория оценивания Теория оценивания Раздел статистики, рассматривающий подходы и методы оценивания некоторых величин или зависимостей, не наблюдаемых в ходе эксперимента непосредственно, по имеющемуся набору данных. Классификация и кластеризация Определение признака принадлежности элемента экспериментального набора данных некоторой (иногда неявно заданной) группе. Задачи являются развитием теории оценивания. Аппроксимация экспериментальных данных Fitting весьма широко употребляется и в теории оценивания и, по сути, означает подгон под какой-то набор имеющихся данных, т.е. подбор функции (через подбор ее параметров) наилучшим образом описывающией характер изменения этого набора данных. 3 / 90 Теория оценивания Что можно оценивать? • стохастические характеристики; (например, E [ x ], VAR [x], COV [x], и др.) • параметры функций распределения; (например для нормального: µx , σx2 ) • функции распределения в графическом виде; (непараметрическое оценивание) • параметры функций, описывающих набор данных; (аппроксимация экспериментальных данных - curve fitting) • сигналы (изменение вектора состояния как функции времени); (фильтр Калмана) В отличии от идентификации систем в теории оценивания речь идет о работе только с некоторым имеющимся набором данных, без какой либо предварительной подготовки эксперимента. 4 / 90 Теория оценивания Оценивание По имеющемуся набору данных численно определить или охарактеризировать статистическую модель. Практически всегда оценивание подразумевает решения некоторой (в зависимости от метода) оптимизационной задачи Статистическая модель Математическая модель, как правило, задаваемая в виде уравнений, описывающая связь одной или нескольких случайных переменных с другими переменными (случайными и детерминистическими). Виды статистических моделей: • непараметрические; • параметрические. Вид модели определяет соответствующий подход оценивания 5 / 90 Непараметрическое оценивание Теория оценивания Непараметрическое оценивание Непараметрическое оценивание Методы оценивания (получение описания одной или нескольких случайных переменных), в которых либо не учитывается закон распределения полученных данных, либо учитывается с неявно определенными параметрами. Другими словами, в методах непараметрического оценивания статистическая модель (и ее структура) не имеет фиксированного числа параметров. 7 / 90 Теория оценивания Непараметрическое оценивание • построение графика функции плотности распределения вероятности; (например, гистограммой) • Kernel Density Estimation; (ядерная оценка (сглаживание) плотности распределения) • непераметрическая регрессия; (на базе ядер, сплайнов, вейвлетов и др.) • ряд методов классификации и кластеризации; (например, kNN, SVM) • проверка статистических гипотез; (теория детектирования, в каком-то смысле) Под ядром понимается некотрая весовая функция 8 / 90 Непараметрическое оценивание Гистограмма Гистограмма Для набора из N случайных переменных {X1 , X2 , . . . , XN } количество ci попаданий Xj в i-ый подинтервал [ai−1 , ai ] исходного интервала [a0 , an ] для i = 1, 2, . . . , n (теоретически, [−∞, +∞]) определяется как: N X ci = {1 : Xj ∈ [ai−1 , ai ]} j=1 Тогда кусочно постоянная фукция h(x), называемая нормализованной гистограммой, оценивается следующим обоазом: ci ci ĥ(x) = = N ∆ai N (ai − ai−1 ) Нормализованная гистограмма есть графическая интерпретация функции плотности распределения вероятности 9 / 90 Непараметрическое оценивание Гистограмма 10 / 90 Непараметрическое оценивание Гистограмма Интегральная гистограмма Для набора из N случайных переменных {X1 , X2 , . . . , XN } количество ci попаданий Xj в i-ый подинтервал [a0 , ai ] исходного интервала [a0 , an ] для i = 1, 2, . . . , n (теоретически, [−∞, +∞]) определяется как: N X ci = {1 : Xj ∈ [a0 , ai ]} j=1 Нормализованная интегральная гистограмма h(x), являющаяся функцией распределения вероятности, оценивается следующим обоазом: ci ci ĥ(x) = = N ∆ai N (ai − ai−1 ) 11 / 90 Непараметрическое оценивание Гистограмма Основная сложность построения гистограммы это выбор оптимального количества разбиений n, для чего существует много формул Правило Стёрджеса: n = 1 + log2 (N) ≈ 1 + 3.322 lg(N) Альтернатива #1: √ n= N Альтернатива #2. Правило Райса: n = 2 N 1/3 Альтернатива #3. Формула Скотта: 3.5 σx n = 1/3 N 12 / 90 Непараметрическое оценивание Kernel Density Estimation Kernel Density Estimation Для набора из N независимых и случайно распределенных случайных переменных {X1 , X2 , . . . , XN } форма (огибающая) функции плотности распределения вероятности определяется как: N N x−X  1 X 1 X j f̂KDE (x) = Kh (x − Xj ) = Kh N Nb b j=1 j=1 где: K(·) - ядро, неотрицательная функция интегрируемая в 1; b > 0 - параметр сглаживания; KDE классический пример проблемы сглаживания функции 13 / 90 Непараметрическое оценивание Сравнение гистограммы и Kernel Density Estimation Красным показаны ядра 14 / 90 Непараметрическое оценивание Виды ядер 15 / 90 Непараметрическое оценивание Kernel Density Estimation Параболическое ядро - оптимальное по среднеквадратичной ошибке (MSE - mean square error) Выбор параметра сглаживания не столь очевиден Оптимальный по средней интегральной квадратичной ошибке: Z    2 MISE(b) = E f̂KDE (x) − f(x) dx Расчет невозможен, потому что неизвестна f(x) Для предположительно нормального закона распределения: !1/5 4 σx5 ≈ 1.06 σx N −1/5 b= 3N 16 / 90 Непараметрическое оценивание Влияние параметра сглаживания b 17 / 90 Непараметрическое оценивание ОЧЕНЬ ВАЖНЫЕ НАБЛЮДЕНИЯ Выбор параметров метода оценивания, как правило, не прост и сводится к поиску некоторой золотой середины между уменьшением смещения (ошибки / bias) и ее разбросом (дисперсией) 18 / 90 Непараметрическое оценивание Проверка статистических гипотез Проверка статистических гипотез Класс задач математической статистики, направленный на численное принятие решения о том, удовлетворяет ли статистическая выборка заданной статистической гипотезе. С практической точки зрения, как правило, речь идет о: • описывается ли случайная переменная из имеющейся выборки заданным законом распределения; • принадлежат ли две случайные величины из имеющейся выборки к одному закону распределения. Одним из наиболее используемых является тест Колмагорова-Смирнова 19 / 90 Параметрическое оценивание Параметрическое оценивание Параметрическое оценивание Определение вектора θ̂, являющегося оценкой вектора параметров θ или интервала, в котором они находятся, статистических моделей M c фиксированными структурой и набором параметров, по имеющемуся набору данных ZN . Истинные параметры Вектор значений θ0 , соответствующий параметрам θ либо заложенным в модель M (в этом случае θ0 известен), либо принимаемый таковым (θ0 неизвестен). Ошибка оценки Вектор значений θ̃, показывающий на сколько ошибчна полученная оценка θ̂. Как правило, θ0 неизвестен, а следовательно, вычисление e напрямую невозможно, поэтому ошибку оценки определяют (оценивают) косвенно. 21 / 90 Параметрическое оценивание Точечная оценка Точечная оценка Вектор θ̂ является точечной оценкой параметров θ статистических моделей M c фиксированными структурой и набором параметров, по имеющемуся набору данных ZN , причем θ̂ и θ0 принадлежат некоторому множеству Θ допустимых значений, что математически формулируется как: θ̂ = argmin L(·) θ∈Θ или θ̂ = argmax L(·) θ∈Θ где: LLMMSE (·) - критерий близости оценки к истинному значению. Иногда используют обозначение θ̂(ZN ) Вычисление точечной оценки не является сложной задачей, сложной задачей является сделать это с достаточной степенью точности 22 / 90 Параметрическое оценивание Практические аспекты задачи параметрического оценивания Критерий оценивания Критерий LLMMSE (·), как положительно определенная функция, показывающий близость получаемой в процессе решения оценки θ̂ к значениям истинных параметров θ0 . Метод оценивания По сути, выбор между решением задачи максимизации или минимизации для заданного критерия LLMMSE (·), а также применение (не обязательно) некоторых механизмов преобразования исходных данных, выбор ограничений. θ̂ = argmin L(·) θ∈Θ или θ̂ = argmax L(·) θ∈Θ Численный метод решения задачи оптимизации Численный метод решения поставленной оптимизационной задачи. 23 / 90 Параметрическое оценивание Количественные характеристики точечных оценок Математическое ожидание оценки Математическое ожидание оценки выражается в математическом ожидании каждого i-го элемента вектора θ̂ по полученному множеству из K оценок θ̂1 , . . . , θ̂K как: h i h i mθ = E θ̂ = E θ̂k (i) k Дисперсия оценки Мера разброса вектора оценок θ̂ по полученному множеству из K оценок θ̂1 , . . . , θ̂K : h i h   i k k Pθ = VAR θ̂ = E θ̂ − mθ θ̂ − mθ > k Pθ - ковариационная матрица Дисперсию можно считать поэлементно 24 / 90 Параметрическое оценивание Неравенство Крамера-Рао Неравенство Крамера-Рао В общем случае неравенство Крамера-Рао показывает нижнюю границу дисперсии оценки θ̂ через информацию Фишера. Другими словами, является средством определения минимально достижимой величины дисперсии оценки, а, как следствие, определения и минимально достижимой точности оценки. Неравенство определено как: 1 Pθ > I(θ) Неравенство справедливо для ряда ограничений на модель Информация Фишера Показывает сколько информации заложено в наборе данных ZN относительно неизвестной величины θ. 25 / 90 Параметрическое оценивание Свойства точечных оценок Смещенность и несмещенность Смещена или нет получаемая оценка θ̂ относительно истинных параметров θ0 на некоторую постоянную величину. Состоятельность Сходится ли оценка θ̂ к истинным параметрам θ0 при объемах данных ZN c N → ∞. Эффективность Является ли полученная оценка θ̂ оптимальной, в контексте минимальной дисперсии, при любых объемах данных ZN . Асимптотическая эффективность Является ли полученная оценка θ̂ оптимальной, в контексте минимальной дисперсии, только при объемах данных ZN c N → ∞. 26 / 90 Параметрическое оценивание Смещенность и несмещенность оценки Смещенность и несмещенность Оценка θ̂ называется несмещенной, если ее математическое ожидание равно истинным значениям θ0 , т.е.: h i E θ̂ = θ0 В противном случае, оценка считается смещенной. На практике речь идет скорее о характеристике метода. Дисперсия несмещенной оценки примет вид: h h i   i k k Pθ = VAR θ̂ = E θ̂ − θ0 θ̂ − θ0 > k что, по сути, означает дисперсию непосредственно ошибки оценки, а Pθ становится ковариационной матрицей ошибок оценок. Для несмещенной оценки можно проще оценить ошибку оценки 27 / 90 Параметрическое оценивание Пример смещенной оценки Выборочная дисперсия: S2N = N 2 1 X  Xi − E [ X ] N i=1 Оценка выборочной дисперсии: " # N 2  2 1 X  E SN = E Xi − E [ X ] N i=1 " # N

2 1 X  =E Xi − mX − E [ X ] − mX N i=1 h
2 i E [ X ] − mX = σX2 − E = σX2 − N−1 2 1 2 σX = σX < σX2 N N 28 / 90 Параметрическое оценивание Пример несмещенной оценки Выборочная исправленная дисперсия: N  2 X 1 S2N = Xi − E [ X ] N−1 i=1 Оценка выборочной исправленной дисперсии: # " N  2 X  2 1 E S =E Xi − E [ X ] N−1 i=1 " # N  X

2 1 =E Xi − mX − E [ X ] − mX N−1 i=1   h
2
2 i N 1 N = E Xi − mX − E E [ X ] − mX N−1 N N−1 N N 1 = σX2 − · σX2 = σX2 N−1 N−1 N 29 / 90 Параметрическое оценивание Эффективность оценки Эффективность оценки Оценка θ̂ называется эффективной, если она является несмещенной и ее дисперсия минимальна (совпадает с нижней границей неравенства Крамера-Рао) среди других оценок, т.е.: h i 1 VAR θ̂ = Pθ = I(θ) Важным следствием этого является то, что ковариационная матрица оценок (ошибок оценок, т.к. она несмещенная) является обратной к информационной матрице Фишера. Асимптотическая эффективность оценки При небольших объемах набора данных ZN оценка θ̂ может быть неэффективной, но при увеличении объема данных ZN , т.е. при N → ∞, θ̂ будет обладать свойствами эффективной оценки, т.е. может считаться асимтотической эффективностой оценкой. 30 / 90 Параметрическое оценивание Состоятельность оценки Состоятельность оценки Оценка θ̂ называется состоятельной, если она сходится к истинным значениям θ0 при объеме набора данных ZN , стремящимся к бесконечности, т.е. при N → ∞. Признаки состоятельной оценки: • сходится к θ0 в среднеквадратичном смысле; • асимптотически эффективная оценка с дисперсией → 0. 31 / 90 Параметрическое оценивание Интервальная оценка Интервальная оценка Две оценки θ̂1 и θ̂2 , полученные по имеющемуся набору данных ZN , называют интервальной оценкой параметров θ статистической модели M с надежностью α, если: h i P θ̂1 6 θ̂ 6 θ̂2 = 1 − α
В этом случае интервал θ̂1 , θ̂2 называется доверительным интервалом, а величина 1 − α доверительной вероятностью (надежностью). Существуют и другие способы получения интервальной оценки, т.е. другие способы формирования интервала 32 / 90 Параметрическое оценивание Важность интервальной оценки На практике вычисляют K оценок θ̂1 , . . . , θ̂K (по наборам данных ZN , полученных экспериментально и с применением методов моделирования Монте-Карло), по которым и определяют доверительный интервал Если дисперсия точечной оценки говорит о том, как точно мы вообще нашли оценку, то интервальная оценка, на практике, относится к тому, как точно мы определили (каков разброс) каждого элемента вектора θ Последнее особенно актуально для динамических систем, где, по факту, после идентификации мы должны учитывать возможный разброс параметров, т.е. рассматривать систему с распределенными параметрами. Подобная интервальная оценка как раз и позволяет сформировать диапазон изменений значений каждого параметра системы 33 / 90 Параметрическое оценивание Общие требования к параметрическому оцениванию Оценка должна: • иметь нулевое смещение (bias); • иметь минимульную дисперсию (variance); • иметь минимульную ошибку оценки; При несмещенной оценки, минимизация ошибки оценки, по сути, означает и минимизацию дисперсии самой оценки Одним из основных подходов является формирования критерия оценивания как величины квадратичной ошибки оценки 34 / 90 Критерий среднеквадратичной ошибки Критерий MSE Критерий среднеквадратичной ошибки (Mean Squared Error) может быть сформулирован для скалярного случая как: h i    h 2
2 i LMSE θ̂ = E θ̂ − θ0 = VAR θ̂ + Bias θ̂, θ0 а для векторного случая как:    h i   LMSE θ̂ = tr VAR θ̂ + Bias θ̂, θ0 2 Минимизация среднеквадратичной ошибки тоже самое, что и минимизация дисперсии (при несмещенной оценке) 35 / 90 Критерий среднеквадратичной ошибки Пример Оценка средней величины: m̂X = N 1 X Xi = m0 N i=1 MSE средней величины:   h
2 i LMSE m̂X = E m̂X − m0 = σ √X N !2 = σX2 N Для распределения Гаусса такая оценка становится лучшей несмещенной оценкой (best unbiased estimate) Метод с несмещенной оценкой, в котором дисперсия оценки минимальна называется минимально-дисперсным несмещенным методом оценивания (minimum-variance unbiased estimator) 36 / 90 Критерий среднеквадратичной ошибки Пример Оценка смещенной величины дисперсии:  2 N 1 X 2 S = Xi − mX N i=1 MSE смещенной величины дисперсии: h  
2 i 2 N − 1 4 = S2 − σX2 LMSE S2 = E σX N2 Оценка несмещенной величины дисперсии:  2 N 1 X 2 SN = Xi − mX N−1 i=1 MSE несмещенной величины дисперсии:   h
2 i LMSE S2N = E S2N − σX2 = 2 σ4 N−1 X 37 / 90 Критерий среднеквадратичной ошибки Метод оценивания путем минимизации среднеквадратичной ошибки Метод MMSE Метод оценивания, в котором минимизируется критерий MSE называется Minimum Meas Squared Error - MMSE, т.е. метод оценивания путем минимизации среднеквадратичной ошибки. Свойства метода оценивания MMSE 1. при конечных значениях среднего и дисперсии   θ̂MMSE = E θ|ZN 2. несмещенная при указанных выше допущениях h i    = E[θ] E θ̂MMSE = E E θ|ZN 3. асимптотически несмещенная  √ 
d N θ̂MMSE − θ0 → N 0, I −1 (θ) 38 / 90 Критерий среднеквадратичной ошибки Свойства оценок полученных по методу MMSE 4. метод оценивания, где оценка сформирована в форме θ̂ = g(ZN ) оптимален, т.е. соответствует θ̂MMSE = gMMSE (ZN ), тогда и только тогда, когда    N E θ̂MMSE − θ0 g(Z ) = 0 что эквивалентно:  E  θ̂MMSE − θ0  θ̂ =0 и для общего вектороного случая    > E θ̂MMSE − θ0 θ̂ =0 Рассмотренное свойство - свойство ортогональности 39 / 90 Критерий среднеквадратичной ошибки Линейный метод оценивания путем минимизации среднеквадратичной ошибки Метод LMMSE Метод оценивания, в котором θ и y = ZN - совместно гауссовские случайные переменные, а критерий LLMMSE (·) сформирован как MSE, где оценка задана в следующей форме: θ̂LMMSE = W y + b т.е. min LLMMSE (·) при θ̂LMMSE = W y + b W,b Подобная форма задания θ̂ возволяет уйти от сложных  вычислений условного математического ожидания E θ|ZN Нужно получить критерий LLMMSE (·) и алгоритм вычисления W и b 40 / 90 Критерий среднеквадратичной ошибки Линейный метод оценивания путем минимизации среднеквадратичной ошибки По требованию несмещенной оценки: h i E θ̂LMMSE = E [ θ ] E [ W y + b ] = mθ W my + b = mθ Откуда получим оптимальное значение b: b = mθ − W my Оценка с оптимальным значением b: θ̂LMMSE = W y + mθ − W my = W Ошибка оценки в этом случае примет вид:
θ̂LMMSE − θ0 = W y − my + mθ − θ0 = W y − my
y − my + mθ
− θ 0 − mθ
41 / 90 Критерий среднеквадратичной ошибки Линейный метод оценивания путем минимизации среднеквадратичной ошибки По принципу ортогональности: h

i θ̂ − θ0 y − my > 0=E h


>i 0=E W y − my − θ0 − mθ y − my 

> 

 0=W E y − my y − my −E θ 0 − mθ y − my > | {z } | {z } Cy Cθy 0 = W Cy − Cθy Следовательно: W = Cθy Cy−1 , а W > = Cy−1 Cyθ 42 / 90 Критерий среднеквадратичной ошибки Линейный метод оценивания путем минимизации среднеквадратичной ошибки И оценка примет вид: θ̂LMMSE = W y + b = Cθy Cy−1 y − my
+ mθ Т.к. θ̂ случайная величина, то ее ковариация: h

>i θ̂ − mθ θ̂ − mθ Cθ̂ = E 

 =W E y − my y − my > W > | {z } Cy = W Cy W > = Cθy Cy−1 Cy Cy−1 Cyθ = Cθy Cy−1 Cyθ 43 / 90 Критерий среднеквадратичной ошибки Линейный метод оценивания путем минимизации среднеквадратичной ошибки ошибки оценки θ̂ − θ0 :

i θ̂ − θ0 θ̂ − θ0 >



>i θ̂ − θ0 W y − my − θ0 − mθ h

i

i θ̂ − θ0 θ̂ − mθ > θ̂ − θ0 y − my > W > − E {z } = 0 по принципу ортогональности h


>
i = −E W y − my − θ0 − mθ θ̂ − mθ > h 



i =E θ 0 − mθ θ0 − mθ > −W E y − my θ̂ − mθ > | {z } | {z } Ковариация h Ce = E h =E h =E | Cθ Cyθ = Cθ − W Cyθ = Cθ − Cθy Cy−1 Cyθ = Cθ − Cθ̂ | {z } C θ̂ 44 / 90 Критерий среднеквадратичной ошибки Линейный метод оценивания путем минимизации среднеквадратичной ошибки Таким образом получаем:   LLMMSE θ̂ = tr (Ce ) Для вычисления W можно использовать: • QR разложение; • разложение Холецкого; • и др. Если y стационарный в широком смысле стохастический процесс, то метод оценивания становится фильтром Винера-Колмогорова 45 / 90 Критерий среднеквадратичной ошибки LMMSE линейного процесса наблюдения LMMSE линейного процесса наблюдения Определяем не величину оценки зависящую от некоторых измерений, а формируем модель измерений, как функцию оценки. Пусть y = A θ + ω(t) E[ω] = 0 Cθω = 0 тогда E [ y ] = A mθ Cy = A Cθ A> + Cω Cθy = Cθ A> 46 / 90 Критерий среднеквадратичной ошибки LMMSE линейного процесса наблюдения Матрицу W определим как: W = Cθ A>  A Cθ A> + Cω −1 Оценку θ̂ получим в виде:  −1   y − A mθ + mθ θ̂ = Cθ A> A Cθ A> + Cω Ковариацию ошибок оценок Ce получим в виде: Ce = Cθ̂ − Cθ  −1 A Cθ = Cθ − Cθ A> A Cθ A> + Cω Не требуется, чтобы число наблюдений N (размерность вектора y) было больше числа неизвестных (размерность вектора θ) 47 / 90 Критерий среднеквадратичной ошибки LMMSE линейного процесса наблюдения. Альтернативная форма Учитывая что: Cθ A>  A Cθ A> + Cω =  A> Cω−1 A + Cθ−1 −1 Ce =  A> Cω−1 A + Cθ−1 −1 −1 A> Cω−1 Матрицу W определим как:  −1 A> Cω−1 = Ce A> Cω−1 W = A> Cω−1 A + Cθ−1 Оценку θ̂ получим в виде: θ̂ = Ce A> Cω−1  y − A mθ  + mθ Подобная оценка сравнима с взвешенным МНК и оценкой Гаусса-Маркова (BLUE - best linear unbiased estimate) 48 / 90 Критерий среднеквадратичной ошибки LMMSE линейного процесса наблюдения. Связь с обычным МНК LMMSE как обычный МНК Если составляющие ω(t) некоррелируемы, т.е. Cω = σ 2 I получаем W=  A> A −1 A> Обычный метод наименьших квадратов 49 / 90 Критерий среднеквадратичной ошибки LMMSE линейного процесса наблюдения. Связь со взвешенным МНК LMMSE как взвешенный МНК Если дисперсия априорной информации о θ бесконечна, т.е. Cθ−1 = 0 получаем W=  A> Cω−1 A −1 A> Cω−1 Взвешенный метод наименьших квадратов с весовой матрицей Cω−1 50 / 90 Критерий среднеквадратичной ошибки Последовательный LMMSE Оценка в режиме реального времени: • перерасчитывать оценку на каждом шаге, но тогда теряется информация с предыдущих наблюдений; • обновлять значение оценки, полученной на предыдущем шаге, при поступлении новых данных (новых измерений). Последовательный LMMSE Метод оценивания, при котором значение оценки в режиме реального времени обновляется при появлении новых данных. Важным допущением является то, что статистические характеристики θ(t) не меняются с течением времени, т.е. θ(t) - стационарный случайный процесс. 51 / 90 Критерий среднеквадратичной ошибки Последовательный LMMSE Пусть на предыдущем шаге была получена оценка θ̂1 с Ce1 Тогда оценка измерений на текущем шаге: ŷ = A θ̂1 Ошибка оценки измерений на текущем шаге: ỹ = y − ŷ   = A θ − θ̂1 + ω(t) = A e1 + ω(t) Оценка на текущем шаге: θ̂2 = θ̂1 + Cθỹ Cỹ−1 ỹ 52 / 90 Критерий среднеквадратичной ошибки Последовательный LMMSE Учитывая, что: E [ ỹ ] = 0 и θ = θ̂1 + e1 Получим: Cỹ = A Ce1 A> + Cω " #    > Cθỹ = E θ̂1 + e1 − mθ A e1 + ω = Ce1 A> Перепишем оценку на текущем шаге:  −1   y − A θ̂1 θ̂2 = θ̂1 + Ce1 A> A Ce1 A> + Cω Ковариация ошибки оценки на текущем шаге:  −1 Ce2 = Ce1 − Ce1 A> A Ce1 A> + Cω A Ce1 53 / 90 Критерий среднеквадратичной ошибки Последовательный LMMSE Последовательный LMMSE является итерационным процессом, на каждой итерации которого выполняется следующие вычисления Kt+1 = Ce1 A> θ̂t+1 Cet+1  A Ce1 A> + Cω   = θ̂t + Kt+1 y − A θ̂t   = I − Kt+1 A Cet −1 Обобщение последовательного LMMSE на случай нестационарности θ(t) ведет к формированию Фильтра Калмана 54 / 90 Критерий среднеквадратичной ошибки Метод наименьших квадратов - Ordinary Least Squares Метод наименьших квадратов Метод оценки параметров функции f(θ, t), в котором функция потерь L(·) представляет собой интегральную меру квадрата отклонения (ошибки) e функции f(θ, x) от заданной y. N N   X X 2 e = L(θ) = f(θ, ti ) − y(ti ) 2 i=1 i=1 а решение задачи оценивания формулируется как: N X θ̂ = argmin e2 θ∈Θ i=1 или θ̂ = argmin kf(θ, t) − y(t)k2 θ∈Θ 55 / 90 Критерий среднеквадратичной ошибки Линейный метод наименьших квадратов Линейный метод наименьших квадратов Метод оценки параметров линейного процесса наблюдения, т.е.: y = Xθ а функция потерь: L(θ) = m  X i=1 yi − n X Xij θj  2 = ky − X θk2 j=1 где: X - матрица размерностью m × n; θ - вектор оцениваемых параметров размерности n; y - вектор наблюдений размерности m. 56 / 90 Критерий среднеквадратичной ошибки Линейный метод наименьших квадратов Решение задачи: X> X θ̂ = X> y Примет вид: θ̂ =  > X X −1 X> y Подобное справедливо при n линейно независимых столбцах X, иначе задача либо переопределенная, либо недоопределенная 57 / 90 Критерий среднеквадратичной ошибки Взвешенный метод наименьших квадратов - Weighted LS Взвешенный метод наименьших квадратов Метод оценки параметров функции f(θ, t), в котором в функции потерь L(·) добавляются весовые коэффициенты W. N N   X X 2 L(θ) = Wii e = Wii f(θ, ti ) − y(ti ) 2 i=1 i=1 где один из вариантов определения весов: Wii = 1/σi2 и для линейного процесса наблюдения решение задачи: X> W X θ̂ = X> W y примет вид: θ̂ =  X> W X −1 X> y 58 / 90 Метод максимального подобия MLE - Maximum Likelihood Estimation Метод максимального подобия Метод оценивания путем максимизации функции подобия: θ̂ = argmax L(θ) θ∈Θ где в качестве функции подобия (в противовес функции потерь) часто используют логарифмическое подобие:   `(θ) = ln L(θ) или среднее логарифмического подобия h  i ˆ = E ln L(θ) `(θ) 59 / 90 Метод максимального подобия Функции подобия Функция подобия для дискретного случая: L(θ|y) = pθ (y) = P [Y = y] Функция подобия для непрерывного случая: L(θ|y) = fθ (y) Относительная функция подобия: ( ) L(θ|y) p θ: > 100 L(θ̂|y) p% - область подобия При θ скалярном значении получаем еще один вид интервальной оценки 60 / 90 Метод апостериорного максимума MAP - Maximum a posteriori Probability Метод апостериорного максимума Пусть θ - случайная переменная с априорным распределением g(θ), апостериорное распределение можно определить как: f(y|θ) g(θ) θ 7→ f(θ|y) = R f(y|θ) g(θ) dθ Θ где знаменатель (совместное подобие) всегда положителен и не зависит от θ, следовательно не влияет на оптимизацию, следовательно оценка: f(y|θ) g(θ) θ̂ = argmax f(θ|y) = argmax R = argmax f(y|θ) g(θ) f(y|θ) g(θ) dθ θ∈Θ θ∈Θ θ∈Θ Θ Метод оценивания путем максимизации апостериорной вероятности 61 / 90 Информационные критерии и выбор модели Информационный критерий Акайке AIC - Akaike Information Criteria Информационный критерий Акайке Критерий, показывающий качество модели M (моделей) для имеющегося набора данных ZN с учетом количества параметров K: 
 AIC = 2 K − 2 ln max L(·) что позволяет сравнивать несколько моделей, но не позволяет делать вывод о близости модели к истинной. Скорректированный информационный критерий Акайке Для случая небольшой выборки объема N: 2 K (K + 1) AICc = AIC + N−K−1 Модель с AICmin лучше описывает набор данных 63 / 90 Информационный критерий Байеса BIC - Bayesian information criterion Информационный критерий Байеса Критерий, показывающий качество модели M (моделей) для имеющегося набора данных ZN с учетом количества параметров K и количества точек в наборе данных N: 
 BIC = ln(N) K − 2 ln max L(·) В отличии от информационного критерия Акайке, информационный критерий Байеса сильнее “штрафует” увеличение количества параметров модели. 64 / 90 Выбор модели Выбор модели Ввиду того, что статистическая модель M, как правило, неизвестна, у нас всегда имеется несколько моделей кандидатов, из которых, жалетаельно на предварительном этапе выбрать наилучшую. Рассмотренные информационные критерии как раз для этого и используются, но, т.к. истинная модель нам не известна, то используют относительное подобие моделей:   AICmin − AICi exp 2 По сути требуется нахождение золотой середины между количеством параметров (чем меньше, тем лучше с позиции жесткой структуры модели) и качеством оценивания (чем больше тем лучше) 65 / 90 Фильтры Калмана Фильтры Калмана Kalman Filters Фильтры Калмана Семейство методов оценивания состояния x в моделях подобного вида по получаемым измерениям в режиме реального времени: (
xk+1 = f xk , uk + ωk
yk = h xk + νk Виды фильтра Калмана • линейный фильтр Калмана (KF); (гибридный, фильтр Калмана-Бюси) • расширенный фильтр Калмана (EKF); (дискретный, непрерывно-дискретный, непрерывный) (модификации: итерационный, робастный, инвариантный) • сигма-точечный фильтр Калмана (UKF); 67 / 90 Линейный фильтр Калмана Kalman Filter Линейный фильтр Калмана Это оптимальный наблюдатель по текущему измерению для системы: ( xk+1 = Ak xk + Bk uk + ωk yk = Ck xk + Dk uk + νk Стохастические процессы: ωk ∼ N (0, Qk ) νk ∼ N (0, Rk ) Стохастические характеристики:   E ω ω > = Q ← ковариационная матрица   E ν ν > = R ← ковариационная матрица   E ω ν > = 0 ← взаимная ковариационная матрица 68 / 90 Линейный фильтр Калмана Инициализация начальных условий Обозначения: x̂n|m x̂k|k Pn|m Pk|k - априорная оценка в момент времени n по m < n измерениям; апостериорная оценка в момент времени k; априорная ковариация ошибок; апостериорная ковариация ошибок. Инициализация начальных условий: x̂0 = E [ x0 ] = mx0 h

>i P0 = E x0 − mx0 x0 − mx0 = Px0 Фильтр Калмана состоит из двух шагов на каждом такте: прогноз (априорная оценка) и корректировка (апостериорная оценка) 69 / 90 Линейный фильтр Калмана Основной алгоритм формирования оценок Априорная оценка: x̂k|k−1 = Ak x̂k−1|k−1 + Bk uk Pk|k−1 = Ak Pk−1|k−1 A> k + Qk Апостериорная оценка: zk = yk − Ck x̂k|k−1 Sk = Ck Pk|k−1 C> k + Rk −1 Kk = Pk|k−1 C> k Sk x̂k|k = x̂k|k−1 + Kk zk

> Pk|k = I − Kk Ck Pk|k−1 I − Kk Ck + Kk Rk K> k 70 / 90 Линейный фильтр Калмана Оптимальная матрица усиления Критерий оценивания (среднеквадратичная ошибка): xk − x̂k|k
Минимизируется tr Pk|k , где
Pk|k = I − Kk Ck Pk|k−1 I − Kk Ck
> + Kk Rk C> k Что достигается если:

> ∂ tr Pk|k = −2 Ck Pk|k + 2 Kk Sk = 0 ∂ Kk Откуда находим оптимальную матрицу усиления:
> Kk Sk = Ck Pk|k = Pk|k−1 C> k −1 Kk = Pk|k−1 C> k Sk 71 / 90 Линейный фильтр Калмана Упрощение вычислений Если модель и начальные условия истинны, то:     E xk − x̂k|k = E xk − x̂k|k−1 = 0 E [ zk ] = 0   Pk|k = COV xk − x̂k|k   Pk|k−1 = COV xk − x̂k|k−1 Sk = COV [zk ] Если Kk - оптимален, то можно использовать:
Pk|k = I − Kk Ck Pk|k−1 На практике такая оптимальность фильтра Калмана недостижима 72 / 90 Фильтр Калмана-Бюси Kalman-Bucy Filter Фильтр Калмана-Бюси Это полностью непрерывный аналог фильтра Калмана для системы: ( ẋ(t) = A(t) x(t) + B(t) u(t) + ω(t) y(t) = C(t) x(t) + D(t) u(t) + ν(t) Оценка требует интегрирования двух дифференциальных уравнений: ẋ(t) = A(t) x(t) + B(t) u(t) + K(t) z(t) Ṗ(t) = A(t) P(t) + P(t) A> + Q(t) − K(t) R(t) K> (t) где: z(t) = y(t) − C(t) x̂(t) K(t) = P(t) C> (t) R−1 Q(t) и R(t) - матрицы спектральной плотности мощности Уравнение Ṗ(t) = . . . - частный случай матричного уравнения Риккати 73 / 90 Гибридный фильтр Калмана Hybrid Kalman Filter Гибридный фильтр Калмана Это фильтр Калмана для оценки вектора состояния непрерывной системы по дискретным измерениям, т.е. модель системы: ( ẋ(t) = A(t) x(t) + B(t) u(t) + ω(t) yk = Ck xk + νk где: xk = x(tk ) ω(t) ∼ N (0, Q(t)) νk ∼ N (0, Rk ) Инициализация начальных условий: x̂0|0 = E [ x(t0 ) ] P0|0 = VAR [x(t0 )] 74 / 90 Гибридный фильтр Калмана Основной алгоритм формирования оценок Априорная оценка: ẋ(t) = A(t) x(t) + B(t) u(t) Ṗ(t) = A(t) P(t) + P(t) A> + Q(t) при x̂(tk−1 ) = x̂k−1|k−1 при P(tk−1 ) = Pk−1|k−1 Интегрируя два дифференциальных уравнения получаем: x̂k|k−1 = x̂(tk ) Pk|k−1 = P(tk ) Апостериорная оценка (идентична дискретному ФК):
−1 > C> Kk = Pk|k−1 C> k k Pk|k−1 Ck + Rk
x̂k|k = x̂k|k−1 + Kk yk − Ck x̂k|k−1

> + Kk Rk K> Pk|k = I − Kk Ck Pk|k−1 I − Kk Ck k 75 / 90 Линейный фильтр Калмана Проблемы использования Проблемы использования линейного фильтра Калмана: • только для линейных моделей; • если модель не истинна или начальные условия определены не верно, то оценка может сильно отличаться от истины; • сложно определять матрицы Qk и Rk ; (одно из решений - ALS - autocovariance least squares)   • малость Qk может привести к вычислению Pk|k с λ Pk|k < 0; (решение через декомпозицию, например P = L D L> ) • фильтр Калмана-Бюси физически не реализуем. Частичное решение этих проблем - расширенный фильтра Калмана (для нелинейных моделей) 76 / 90 Расширенный фильтр Калмана EKF - Extended Kalman Filter Расширенный фильтр Калмана Метод оценивания для нелинейных моделей вида: (
xk+1 = f xk , uk + ωk
yk = h xk + νk основной идеей которого является линеаризация модели на каждом шаге, т.е. f можно использовать для оценки x, но нельзя использовать напрямую для определения P. Линеаризация происходит за счет вычисления матрица Якоби (получаем расширенный фильтр Калмана первого порядка) 77 / 90 Расширенный фильтр Калмана Алгоритм формирования оценок Априорная оценка: x̂k|k−1 = f xk−1|k−1 , uk Fk = ∂f ∂x
x̂k−1|k−1 , uk Pk|k−1 = Fk Pk−1|k−1 F> k + Qk В случае итерационной линеаризации получаем итерационный расширенный фильтр Калмана (Iterated EKF) EKF более высоких порядков требует вычисление частных производных более высокого порядка 78 / 90 Расширенный фильтр Калмана Алгоритм формирования оценок Апостериорная оценка:
zk = yk − h xk|k−1 Hk = ∂h ∂x x̂k|k−1 Sk = Hk Pk|k−1 H> k + Rk −1 Kk = Pk|k−1 H> k Sk x̂k|k = x̂k|k−1 + Kk zk
Pk|k = I − Kk Hk Pk|k−1 I − Kk Hk
> + Kk Rk K> k EKF не оптимален по определению (оптимален только в случае истинной линейной модели) 79 / 90 Расширенный фильтр Калмана Виды расширенного фильтра Калмана Виды расширенного фильтра Калмана • полностью непрерывный EKF; • гибридный EKF. Все аналогично линейному фильтру Калмана, за исключением нелинейных f(·) и h(·) и линеаризации системы на каждом шаге Для полностью непрерывного EKF: F(t) = ∂f ∂x , H(t) = x̂(t) , u(t) ∂h ∂x x̂(t) Для гибридного EKF: F(t) = ∂f ∂x , x̂(t) , u(t) Hk = ∂h ∂x x̂k|k−1 80 / 90 Сигма-точечный фильтр Калмана UKF - Unscented Kalman Filter Сигма-точечтный фильтр Калмана Метод оценивания для нелинейных моделей вида: (
xk+1 = f xk , uk + ωk
yk = h xk + νk основной идеей которого является линеаризация модели (методом сигма-точечечного преобразования) на каждом шаге. UKF призван повысить устойчивость EKF, но на практике, ввиду линеаризации, все равно возможна неустойчивость алгоритма 81 / 90 Сигма-точечный фильтр Калмана Весовые коэффициенты Вычисление весовых коэффициентов: λ = α2 (L + κ) − L λ W0x = L+λ λ + (1 − α2 + β) W0C = L+λ λ WiC = Wix = 2 (L + λ) где: λ κ α β - масштабный коэффициент; масштабный коэффициент (обычно принимает равным нулю); коэффициент, характеризует разброс сигма-точек относительно x; коэффициент, характеризующий информацию о распределении x; Обычно: β = 2 для нормального распределения, α  1 82 / 90 Сигма-точечный фильтр Калмана Расширение вектора состояния и ковариационой матрицы Расширенный вектор состояния и ковариационная матрица:     xk Pk 0 0     sk = ωk  , Ck =  0 Qk 0  0 0 Rk νk Размерность расширенного вектора состояния - L Инициализация начальных условий:    E x0|0   ŝ0|0 =  E [ ω0 ]  E [ ν0 ]  C0|0 = E s0|0 − ŝ0|0
s0|0 − ŝ0|0
>    P0|0 0   =  0 Q0|0 0  R0|0 83 / 90 Сигма-точечный фильтр Калмана Алгоритм формирования оценок Вычисление   Xk−1|k−1   Xk−1|k−1 i   Xk−1|k−1 2 L + 1 сигма-точек: = sk−1|k−1 = sk−1|k−1 + q q
 L + λ Ck−1|k−1 , i = 1 , . . . , L i 
L + λ Ck−1|k−1 , i = L + 1, ... , 2L = sk−1|k−1 + i−L   где · i - i-ый столбец соответствующей матрицы √ Вычисление матрицы · можно эффективно получить через разложение Холецкого, т.е. B = A A∗ , что для вещественных матриц B = A A i В матричной форме: h i p p X = Xk−1|k−1 Xk−1|k−1 + γ Ck−1|k−1 Xk−1|k−1 + γ Ck−1|k−1 γ= √ L + λ, X - квадратная матрица размерностью 2 L + 1 84 / 90 Сигма-точечный фильтр Калмана Алгоритм формирования оценок Априорная оценка сигма-точек:   x x ω Xk|k−1 = F Xk−1|k−1 , uk−1 , Xk−1|k−1   x ν Yk|k−1 = H Xk|k−1 , Xk−1|k−1 F(·) : R2 L+1 → Rn , H(·) : R2 L+1 → Rq Априорные оценки: 2L   X x x̂k|k−1 = Wix Xk|k−1 i i=0 ŷk|k−1 = 2L X Wix  Yk|k−1 i=0 Pk|k−1 = 2L X i=0 WiC   i x Xk|k−1  i − x̂k|k−1  x Xk|k−1 >  i − x̂k|k−1 85 / 90 Сигма-точечный фильтр Калмана Алгоритм формирования оценок Расчет ковариаций: Pyy = 2L X WiC  x Yk|k−1  x Xk|k−1  i=0 Pxy = 2L X WiC  i=0 i i − ŷk|k−1 − x̂k|k−1   > x Yk|k−1  x Yk|k−1  i − ŷk|k−1 > i − ŷk|k−1 Апостериорные оценки: Kk = Pxy  Pyy −1 x̂k|k = x̂k|k−1 + Kk  yk − ŷk|k−1  > Pk|k = Pk|k−1 − Kk Pyy Kk 86 / 90 Сигма-точечный фильтр Калмана Алгоритм формирования оценок Апостериорная оценка расширенного вектора состояния:   xk|k   ŝk|k = E [ ωk ] E [ νk ] Апостериорная оценка расширенной ковариационной матрицы:   Pk|k 0 0   Ck|k =  0 Rk 0  Qk При E [ ωk ] = 0 и E [ νk ] = 0 алгоритм можно упростить, что разделит вычисление сигма-точек на два этапа 87 / 90 Сигма-точечный фильтр Калмана Случай E [ ωk ] = 0 и E [ νk ] = 0 Инициализация начальных условий:   x̂0|0 = E x0|0 
P0|0 = E x0|0 − x̂0|0 x0|0 − x̂0|0
>  Вычисление сигма-точек: h i p p Xk−1|k−1 = xk−1|k−1 xk−1|k−1 + γ Pk−1|k−1 xk−1|k−1 − γ Pk−1|k−1   ∗ Xk|k−1 = F Xk−1|k−1 , uk−1 88 / 90 Сигма-точечный фильтр Калмана Случай E [ ωk ] = 0 и E [ νk ] = 0 Априорные оценки: 2L X ∗ x̂k|k−1 = Wix Xk|k−1 i=0 Pk|k−1 = 2L X WiC  x Xk|k−1  i=0 i − x̂k|k−1  x Xk|k−1 >  i − x̂k|k−1 + Qk Перерасчет сигма-точек и расчет измерений: h √ √ i ∗ ∗ ∗ Xk|k−1 = Xk|k−1 − γ Qk Xk|k−1 + γ Qk Xk|k−1   Yk|k−1 = H Xk|k−1 ŷk|k−1 = 2L X i=0 Wix  Yk|k−1  i 89 / 90 Сигма-точечный фильтр Калмана Случай E [ ωk ] = 0 и E [ νk ] = 0 Расчет ковариационных матриц и апостериорных оценок: >   2L   X x x C + Rk Yk|k−1 − ŷk|k−1 Yk|k−1 − ŷk|k−1 Pyy = Wi Pxy = 2L X WiC  x Xk|k−1 i=0 Kk = Pxy i i i=0  Pyy  i − x̂k|k−1  x Yk|k−1 >  i − ŷk|k−1 −1 x̂k|k = x̂k|k−1 + Kk  yk − ŷk|k−1  > Pk|k = Pk|k−1 − Kk Pyy Kk Не требуется вычислений апостериорного расширенного вектора состояния и расширенной ковариационной матрицы 90 / 90