Фильтр Калмана для минимизации энтропийного значения случайной погрешности с не Гауссовым распределением01.02.2019 17:34

betrvcijf5u2mdhyv2qtalcjpga.png


Введение


На Habr математическое описание работы фильтра Калмана и особенности его применения рассматривались в следующих публикациях [1÷10]. В публикации [2] в простой и доходчивой форме рассмотрен алгоритм работы фильтра Калмана (ФК) в модели «пространства состояний», Следует отметить, что исследование систем контроля и управления во временной области с помощью переменных состояния широко используется в последнее время благодаря простоте проведения анализа [11].

Публикация [8] представляет значительный интерес именно для обучения. Очень эффективен методический приём автора, который начал свою статью с рассмотрения распределения случайной погрешности Гаусса, рассмотрел алгоритм ФК и закончил простой итерационной формулой для подбора коэффициента усиления ФК. Автор ограничился рассмотрением распределения Гаусса мотивируя это тем, что при достаточно больших $n$ (многократных измерений) закон распределения суммы случайных величин стремится к распределению Гаусса.

Теоретически такое утверждение, безусловно, справедливо, однако на практике число измерений в каждой точке диапазона не может быть очень большим. Сам R.E. Kalman получил результаты о минимуме ковариации фильтра на базе ортогональных проекций, без предположения о гауссовости ошибок измерений [12].

Целью настоящей публикации является исследование возможностей фильтра Калмана для минимизации энтропийного значения случайной погрешности с не Гауссовым распределением.
Для оценки эффективности фильтра Калмана при идентификации закона распределения или суперпозицией законов по экспериментальным данным воспользуемся информационная теорией измерений основанной на теории информации К. Шеннона, согласно которой информация, подобно физической величине, может быть измерена и оценена.
Основное достоинство информационного подхода к описанию измерений состоит в том, что размер энтропийного интервала неопределенности может быть найден строго математически для любого закона распределения. Что устраняет исторически сложившийся произвол, неизбежный при волевом назначении различных значений доверительной вероятности.Это особенно важно и в учебном процессе, когда обучающий может наблюдать уменьшение неопределённости измерений при применении фильтрации Калмана на заданной ему числовой выборке[13,14].

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

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

Плотность распределения вероятностей для каждого столбца гистограммы [14] шириной $d$ равна $p_{i}(x_{i})=n_{i}/(n\cdot d) $, отсюда оценка вероятностей энтропии определяется как $H\left (x\right)= \int_{-\infty }^{+\infty}p\left ( x\right)\ln p\left (x\right )dx$ при нахождении энтропии по гистограмме из $m $ столбцов получим соотношение:

$\displaystyle H\left (x\right )= -\sum_{i=1}^{m}\int_{\tilde{x_{i}}-\frac{d}{2}}^{\tilde{x_{i}}+\frac{d}{2}}\frac{n_{i}}{nd}\ln\frac{n_{i}}{nd} =\sum_{i=1}^{m}\frac{n_{i}}{n}\ln \frac{n}{n_{i}}+\ln d$

Представим выражение для оценки энтропии в виде:

$H\left (x\right )= \ln \left [ d\prod_{i=1}^{m}\left (\frac{n}{n_{i}} \right)^{\frac{n_{i}}{n}}\right ]$

Получим выражение для оценки энтропийного значения случайной величины:

$\Delta _{e}= \frac{1}{2}e^{H\left (x\right )}=\frac{dn}{2}10^{-\frac{1}{n}\sum_{1}^{m}n_{i}\lg n_{i}} $

Классификацию законов распределения осуществляют на плоскости в координатах энтропийного коэффициента $k=\frac{\Delta _{e}}{\sigma}$ и контрэксцесса $\psi =\frac{\sigma ^{2}}{\sqrt{\mu _{4}}}$ где $\mu _{4}=\frac{1}{n}\sum_{1}^{n}\left ( x_{i}-\bar{X} \right )^{4}$.

Для всех возможных законов распределения \psi изменяется от 0 до 1, а k от 0 до 2,066, поэтому каждый закон может характеризоваться некоторой точкой. Покажем это, используя следующую программу:

Плоскость законов распределения 
import matplotlib.pyplot as plt
from numpy import*
from scipy.stats import *
def graf(a):#группировка данных
         a.sort()
         n=len(a)
         m= int(10+sqrt(n))
         d=(max(a)-min(a))/m
         x=[];y=[]
         for i in arange(0,m,1):
                  x.append(min(a)+d*i)
                  k=0
                  for j in a:
                           if min(a)+d*i <=j


8m06t21wrfsx_z5xbtvdjwfjhgs.png

На плоскости в координатах $k,\psi$ удалены от остальных распределений, распределения Парето, Коши хотя они и относятся к разным областям применения первое в физике, а второе в экономике. Для сравнения следует выбрать находящееся в вершине классификации нормальное распределение Гаусса. Все приведенные ниже сравнения выполняются на ограниченной выборке и носят характер демонстрации возможностей ФК на примере численного определения энтропийной погрешности.

Выбор алгоритма фильтра Калмана


В каждой выбранной точке диапазона измерений проводятся многократные измерения и их результат сравнивается с мерой которую ФК «не знает». Поэтому следует выбрать ФК, например Kalman Filter to Estimate a Constant [16]. Однако я предпочитаю Python и остановился на варианте [16] с обширной документацией. Приведу описание алгоритма:

Поскольку константа всегда одна модель системы можно представить в виде:

$x_{k}=x_{k-1}+w_{k}$, (1)

Для модели матрица перехода вырождается в единицу, а матрица управления в ноль.Модель измерения примет вид:

$y_{k}=y_{k-1}+v_{k}$, (2)

Для модели (2) матрица измерений превращаются в единицу, а ковариационные матрицы $P,Q,R$ превращаются в дисперсии.На очередном $k$-м шаге, до прихода результатов измерения, скалярный фильтр Калмана пытается по формуле (1) оценить новое состояние системы:

$\hat{x}_{k/(k-1)}=\hat{x}_{(k-1)/(k-1)}$, (3)

Уравнение (3) показывает, что априорная оценка на следующем шаге равна апостериорной оценке, сделанной на предыдущем шаге. Априорная оценка дисперсии ошибки:

$ P_{k/(k-1)}=P_{(k-1)/(k-1)}+Q_{k}$, (4)

По априорной оценке состояния $\hat{x}_{k/(k-1)} $ можно вычислить прогноз измерения:

$\hat{y}_{k}=\hat{x}_{k/(k-1)}$, (5)

После того, как получено очередное измерение $ y_{k}$, фильтр рассчитывает ошибку своего прогноза $k$-го измерения:

$ e_{k}=y_{k}-\hat{y}_{k}=y_{k}-\hat{x}_{k/(k-1)} $, (6)

Фильтр корректирует свою оценку состояния системы, выбирая точку, лежащую где-то между первоначальной оценкой $\hat{x}_{k/(k-1)} $и точкой, соответствующей новому измерению $y_{k}$:

$e_{k}=y_{k}-\hat{y}_{k}=y_{k}-\hat{x}_{k/(k-1)}$, (7)

где $G _{k} $— коэффициент усиления фильтра.

Также корректируется оценка дисперсии ошибки:

$P_{k/(k)}=(1-G_{k})\cdot P_{k/(k-1)}$, (8)

Можно доказать, что дисперсия ошибки $e_{k}$ равна:

$S_{k}=P_{k/(k-1)}+R_{k}$, (9)

Коэффициент усиления фильтра, при котором достигается минимальная ошибка оценки состояния системы определяется из соотношения:

$G_{k}=P_{k/(k-1)}/S_{k}$, (10)

Минимизация энтропийной погрешности ФК для шумов с распределением Коши, Парето и Гаусса

1. В теории вероятности плотность распределения Коши определяется из соотношения:

$f(x)=\frac{b}{\pi\cdot (1-x^{2})}$

Для этого распределения невозможно оценить погрешность методами теории вероятности ($\sigma =\infty$), но информационная теория позволяет это сделать:

Программа минимизации ФК энтропийной погрешности от шумов Коши 
from numpy import *
import matplotlib.pyplot as plt
from scipy.stats import *
def graf(a):#группировка данных для определения энтропийной погрешности
         a.sort()
         n=len(a)
         m= int(10+sqrt(n))
         d=(max(a)-min(a))/m
         x=[];y=[]
         for i in arange(0,m,1):
                  x.append(min(a)+d*i)
                  k=0
                  for j in a:
                           if min(a)+d*i <=j


p9dag8khd79xznepwepwoygowxo.png

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

2. В теории вероятности плотность распределения Парето с параметрами $x_{m}$ и $k$ определяется из соотношения:

$f_{X}(x)= \left\{\begin{matrix} \frac{kx_{k}^{m}}{x^{k+1}},& x\geq x_{m}\\ 0,& x< x_{m} \end{matrix}\right. $

Следует отметить, что распределение Парето встречается не только в экономике. Можно привести следующий пример распределение размера файла в интернет-трафике по TCP-протоколу.

3. В теории вероятности плотность нормального распределения (Гаусса) с математическим ожиданием $\mu$ и среднеквадратическим отклонением $\sigma$ определяется из соотношения:

$f(x)=\frac{1}{\sigma \sqrt{2\pi}}\cdot e^{-\frac{(x-\mu)^{2}}{2\sigma ^{2}}}$

Определение минимизации ФК энтропийной погрешности от шумов с распределением Гаусса приводим для сравнения с не Гауссовыми распределениями Коши и Парето.

Программа минимизации ФК энтропийной погрешности от шумов нормального распределения 
from numpy import *
import matplotlib.pyplot as plt
from scipy.stats import *
def graf(a):#группировка данных для определения энтропийной погрешности
         a.sort()
         n=len(a)
         m= int(10+sqrt(n))
         d=(max(a)-min(a))/m
         x=[];y=[]
         for i in arange(0,m,1):
                  x.append(min(a)+d*i)
                  k=0
                  for j in a:
                           if min(a)+d*i <=j


bb7q_wlwdsbp6cgos62axp4afsa.png

Распределение Гаусса обеспечивает более высокую стабильность результата при 50 измерениях и для приведенного графика энтропийная погрешность уменьшается в 2,2 раза.

Минимизация ФК энтропийной погрешности по выборке экспериментальных данных с неизвестным законом распределения шумов

Программа минимизации ФК энтропийной погрешности ограниченной выборки экспериментальных данных 
from numpy import *
import matplotlib.pyplot as plt
from scipy.stats import *
def graf(a):#группировка данных для определения энтропийной погрешности
         a.sort()
         n=len(a)
         m= int(10+sqrt(n))
         d=(max(a)-min(a))/m
         x=[];y=[]
         for i in arange(0,m,1):
                  x.append(min(a)+d*i)
                  k=0
                  for j in a:
                           if min(a)+d*i <=j


e1dtcui8l759r5tuybvkseulzi0.png

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

Вывод


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

Ссылки

© Habrahabr.ru