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

  • Tutorial

Введение


На 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<min(a)+d*(i+1):
                                    k=k+1
                  y.append(k)
         h=(1+0.5*m/n)*0.5*d*n*10**(-sum([w*log10(w) for w in y if w!=0])/n)     
         k=h/std (a)
         mu4=sum ([(w-mean (a))**4 for w in a])/n
         psi=(std(a))**2/sqrt(mu4)
         return  psi,k
b=800#объём контрольной выборки
plt.title('Плоскость законов распределения', size=12)
plt.xlabel('Коэффициент $\psi$', size=12)
plt.ylabel('Энтропийный коэффициент k', size=12)
a=uniform.rvs( size=b)
psi,k,=graf(a)
nr="Равномерное : k=%s,$\psi$=%s "%(round(k,3),round(psi,3))
plt.plot(psi,k,'o',  label=nr)
a=logistic.rvs( size=b)
psi,k,=graf(a)
nr="Логистическое:k=%s,$\psi$=%s "%(round(k,3),round(psi,3))
plt.plot(psi,k,'o',  label=nr)
a=norm.rvs( size=b)
psi,k,=graf(a)
nr="Нормальное :K=%s,$\psi$i=%s "%(round(k,3),round(psi,3))
plt.plot(psi,k,'o',   label=nr)
a= erlang.rvs(4,size=b)
psi,k,=graf(a)
nr=" Эрланга :k=%s,$\psi$=%s "%(round(k,3),round(psi,3))
plt.plot(psi,k, 'o',  label=nr)
a= pareto.rvs(4,size=b)
psi,k,=graf(a)
nr="Парето :k=%s,$\psi$=%s "%(round(k,3),round(psi,3))
plt.plot(psi,k, 'o',  label=nr)
a = cauchy.rvs(size=b)
psi,k,=graf(a)
nr="Коши :k=%s,$\psi$=%s "%(round(k,3),round(psi,3))
plt.plot(psi,k,'o',   label=nr)
c = 0.412
a = genlogistic.rvs(c, size=b)
psi,k,=graf(a)
nr="Логистическое -1:k=%s,$\psi$=%s "%(round(k,3),round(psi,3))
plt.plot(psi,k,'o',   label=nr)
mu=0.6
a = poisson.rvs(mu, size=b)
psi,k,=graf(a)
nr="Пуассона:k=%s,$\psi$=%s "%(round(k,3),round(psi,3))
plt.plot(psi,k,'o',   label=nr)
a= laplace.rvs(size=b)
psi,k,=graf(a)
nr="Лапласа:k=%s,$\psi$=%s "%(round(k,3),round(psi,3))
plt.plot(psi,k,'o',   label=nr)
plt.legend(loc='best')
plt.grid()
plt.show()   




На плоскости в координатах $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<min(a)+d*(i+1):
                                    k=k+1
                  y.append(k)
         h=(1+0.5*m/n)*0.5*d*n*10**(-sum([w*log10(w) for w in y if w!=0])/n)         
         return  h
def Kalman(a,x,sz,R1):
    R = R1*R1  # Дисперсия    
    Q = 1e-5 # Дисперсия случайной величины в модели системы 
    # Выделение памяти под массивы:
    xest1 = zeros(sz)  # Априорная оценка состояния
    xest2 = zeros(sz)  # Апостериорная оценка состояния  
    P1 = zeros(sz)     # Априорная  оценка ошибки  
    P2 = zeros(sz)     # Апостериорная оценка ошибки  
    G = zeros(sz)      # Коэффициент усиления фильтра
    xest2[0] = 0.0
    P2[0] = 1.0
    for k in arange(1, sz): # Цикл по отсчётам времени. 
             xest1[k] = xest2[k-1] # Априорная оценка состояния. 
             P1[k] = P2[k-1] + Q# Априорная оценка ошибки.
             # После получения нового значения измерения вычисляем апостериорные оценки : 
             G[k] = P1[k] / ( P1[k] + R )
             xest2[k] = xest1[k] + G[k] * ( a[k] - xest1[k] )
             P2[k] = (1 - G[k]) * P1[k]
    return  xest2,P1
nr="распределением Коши "
x =2# Истинное значение измеряемой величины фильтру неизвестно)
R1 = 0.1  # Ср . кв . ошибка измерения . 
sz = 50   # Число итераций.
a = cauchy.rvs(x, R1, size=sz)
xest2,P1=Kalman(a,x,sz,R1)
plt.plot(a, 'k+', label='Зашумлённые измерения')
plt.plot(xest2,'b-', label='Апостериорная оценка')
plt.axhline(x, color='g', label='Истинное значение')
plt.axis([0, sz,-x, 2*x])
plt.legend()
plt.xlabel('Номер итерации')
plt.ylabel(u'Измеряемая величина')
h1=graf(a)
h2=graf(xest2)
plt.title('Подавление шумов с %s.\n Энтропийная погрешность до ФК $\Delta $1=%s после ФК $\Delta $2=%s '%(nr,round(h1,3),round(h2,3)), size=12)
plt.grid()
plt.show()




Вид графика может изменяться как при перезагрузке программы(новой генерации выборки распределения), так и в зависимости от числа измерений и параметров распределения, но одно останется неизменным ФК минимизирует значение энтропийной погрешности вычисленной на основе информационной теории измерений. Для приведенного графика ФК снижает энтропийную погрешность в 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<min(a)+d*(i+1):
                                    k=k+1
                  y.append(k)
         h=(1+0.5*m/n)*0.5*d*n*10**(-sum([w*log10(w) for w in y if w!=0])/n)         
         return  h
def Kalman(a,x,sz,R1):
    R = R1*R1  # Дисперсия    
    Q = 1e-5 # Дисперсия случайной величины в модели системы 
    # Выделение памяти под массивы:
    xest1 = zeros(sz)  # Априорная оценка состояния
    xest2 = zeros(sz)  # Апостериорная оценка состояния  
    P1 = zeros(sz)     # Априорная  оценка ошибки  
    P2 = zeros(sz)     # Апостериорная оценка ошибки  
    G = zeros(sz)      # Коэффициент усиления фильтра
    xest2[0] = 0.0
    P2[0] = 1.0
    for k in arange(1, sz): # Цикл по отсчётам времени. 
             xest1[k] = xest2[k-1] # Априорная оценка состояния. 
             P1[k] = P2[k-1] + Q# Априорная оценка ошибки.
             # После получения нового значения измерения вычисляем апостериорные оценки : 
             G[k] = P1[k] / ( P1[k] + R )
             xest2[k] = xest1[k] + G[k] * ( a[k] - xest1[k] )
             P2[k] = (1 - G[k]) * P1[k]
    return  xest2,P1
nr="распределением Гаусса "
x =2# Истинное значение измеряемой величины фильтру неизвестно)
R1 = 0.1  # Ср . кв . ошибка измерения . 
sz = 50   # Число итераций.
a=norm.rvs( x, R1, size=sz)
xest2,P1=Kalman(a,x,sz,R1)
plt.plot(a, 'k+', label='Зашумлённые измерения')
plt.plot(xest2,'b-', label='Апостериорная оценка')
plt.axhline(x, color='g', label='Истинное значение')
plt.axis([0, sz,-x, 2*x])
plt.legend()
plt.xlabel('Номер итерации')
plt.ylabel(u'Измеряемая величина')
h1=graf(a)
h2=graf(xest2)
plt.title('Подавление шумов с %s.\n Энтропийная погрешность до ФК $\Delta $1=%s после ФК $\Delta $2=%s '%(nr,round(h1,3),round(h2,3)), size=12)
plt.grid()
plt.show()




Распределение Гаусса обеспечивает более высокую стабильность результата при 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<min(a)+d*(i+1):
                                    k=k+1
                  y.append(k)
         h=(1+0.5*m/n)*0.5*d*n*10**(-sum([w*log10(w) for w in y if w!=0])/n)
         k=h/std (a)
         mu4=sum ([(w-mean (a))**4 for w in a])/n
         psi=(std(a))**2/sqrt(mu4)
         return  h
def Kalman(a,x,sz,R1):
    R = R1*R1  # Дисперсия    
    Q = 1e-5 # Дисперсия случайной величины в модели системы 
    # Выделение памяти под массивы:
    xest1 = zeros(sz)  # Априорная оценка состояния
    xest2 = zeros(sz)  # Апостериорная оценка состояния  
    P1 = zeros(sz)     # Априорная  оценка ошибки  
    P2 = zeros(sz)     # Апостериорная оценка ошибки  
    G = zeros(sz)      # Коэффициент усиления фильтра
    xest2[0] = 0.0
    P2[0] = 1.0
    for k in arange(1, sz): # Цикл по отсчётам времени. 
             xest1[k] = xest2[k-1] # Априорная оценка состояния. 
             P1[k] = P2[k-1] + Q# Априорная оценка ошибки.
             # После получения нового значения измерения вычисляем апостериорные оценки : 
             G[k] = P1[k] / ( P1[k] + R )
             xest2[k] = xest1[k] + G[k] * ( a[k] - xest1[k] )
             P2[k] = (1 - G[k]) * P1[k]
    return  xest2,P1

R1 = 0.9  # Ср . кв . ошибка измерения . 
a=[ 0.203, 0.154, 0.172, 0.192, 0.233, 0.181, 0.219, 0.153, 0.168, 0.132, 0.204, 0.165, 0.197, 0.205, 0.143,   0.201, 0.168, 0.147, 0.208, 0.195, 0.153, 0.193, 0.178, 0.162, 0.157, 0.228, 0.219, 0.125, 0.101, 0.211,0.183, 0.147,    0.145, 0.181,0.184, 0.139, 0.198, 0.185, 0.202, 0.238, 0.167, 0.204, 0.195, 0.172, 0.196, 0.178, 0.213, 0.175, 0.194,   0.178, 0.135, 0.178, 0.118, 0.186, 0.191]
sz = len(a)   # Число итераций
x =0.179# Истинное значение измеряемой величины фильтру неизвестно)
xest2,P1=Kalman(a,x,sz,R1)
plt.plot(a, 'k+', label='Зашумлённые измерения')
plt.plot(xest2,'b-', label='Апостериорная оценка')
plt.axhline(x, color='g', label='Истинное значение')
plt.axis([0, sz,-x, 2*x])
plt.legend()
plt.xlabel('Номер итерации')
plt.ylabel('Измеряемая величина')
h1=graf(a)
nr=" неизвестным распределением"
h2=graf(xest2)
plt.title('Подавление шумов с %s \n Энтропийная погрешность до ФК $\Delta $1=%s после ФК $\Delta $2=%s '%(nr,round(h1,3),round(h2,3)), size=12)
plt.grid() 
plt.show()




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

Вывод


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

Ссылки
1. Неортогональная БИНС для малых БПЛА
2. Фильтр Калмана — !cложно?
3. Без запаха фильтрации и нелинейного оценивания *
4. На пороге дополненной реальности: к чему готовиться разработчикам (часть 2 из 3)
5. Классическая механика: о диффурах «на пальцах»
6. Фильтр Калмана — Введение
7. Генератор Федеративного Фильтра Калмана с использованием Генетических Алгоритмов.
8. Фильтр Калмана
9. Использование фильтра Калмана для определения производных измеряемой величины
10. Простая модель адаптивного фильтра Калмана средствами Python
11. Пространство состояний в задачах проектирования систем оптимального управления
12. A New Approach to Linear Filtering and Prediction Problems1
13. Вероятностный и информационный анализ результатов измерений на Python
14. Подбор закона распределения случайной величины по данным статистической выборки средствами Python
15. Kalman filtering
16. Kalman Filter to Estimate a Constant
Поделиться публикацией

Комментарии 17

    0
    1. Не могли бы вы детально рассказать как вы перешли от дискретной энтропии H=-sum(p*ln(p)) с суммой безразмерных величин, к интегрированию плотности распределения случайной величины с размерными величинами под логарифмом?
    2. Почему у вас постоянно пропадают знаки, то энтропии, то в распределении Коши.
    3. И какое отношение «энтропийная погрешность» имеет отношение к погрешности измерения?
    4. Не кажется ли вам что в учебных материалах, надо определять все величины по ходу пьесы, а не ссылать в список литературы без приведения минимально информации из неё необходимой для дальнейшего объяснения.
    5. Что бы «учебную задачу данной статьи можно считать выполненной» эту задачу следовало бы хотя бы определить. А так получили какие то не ведомые цифры ни имеющие вообще никакого отношения ни к какой практической задаче.

    ps: Ваше рвение писать наукообразные статьи похвально. Но пока получается не очень.
      –1
      «1. Не могли бы вы детально рассказать как вы перешли от дискретной энтропии H=-sum(p*ln(p)) с суммой безразмерных величин, к интегрированию плотности распределения случайной величины с размерными величинами под логарифмом?»

      Рассказывать о том что относительные величины не имеют размерности не нужно, это и так понятно.!

      «2. Почему у вас постоянно пропадают знаки, то энтропии, то в распределении Коши.»

      Правило логарифма отношения величин изучают в школе ln(a/b)=-ln(b/a)

      «3.И какое отношение «энтропийная погрешность» имеет отношение к погрешности измерения?»

      Энтропийная погрешность измерений одна из форм представления погрешности измерений. Школьная программа по информатике.

      «4.Не кажется ли вам что в учебных материалах, надо определять все величины по ходу пьесы, а не ссылать в список литературы без приведения минимально информации из неё необходимой для дальнейшего объяснения.»

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

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

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

      Снижение погрешность измерений самая актуальная практическая задача. При этом какие то неведомые цифры( по Вашему определению)" характеризуют энтропийную погрешность в полном соответствии с информационной теорией измерений.

      ps. Ваше рвение писать наукообразные комментарии похвально. Но пока получается не очень.
        0
        1. Плотность распределения имеет размерность 1/dx
        2. К логарифмам претензий нет.
        image
        image
        3. Энтропия это мера количества информации/состояний системы. А что такое энтропийная погрешность в вашем случае.
        4. Если вы пишете учебный материал. Те кто его пытаются читать должны получать информацию последовательно, без разрывов. (Хорошо бы еще полезную). Так что бы не нужно было лазить по ссылкам что бы понять что обозначают буквы которыми вы оперируете. Из ссылок надо приводить минимальную информацию, достаточную для дальнейшего чтения вашего материала без лазания по ссылке.
        «Методика расчета энтропийной погрешности приведена» — где она приведена, только в коде на питоне где вы разбиваете область значений на 17 интервалов?
        Интересует не как определена энтропийная погрешность (математически не корректно), а какой смысл в её вычислении и почему вы её пытаетесь использовать.

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

        Извините какая цель? Убрать плац с помощью лома? Тут важен сам процесс?

        Снижение погрешность измерений самая актуальная практическая задача.

        Вы можете оценить погрешность измерения. Оценить распределение ошибок и т.п.
        Но практическая задача снижения погрешности лежит совсем в другой области. Для более точного измерения величины нужны более точные приборы, большая статистика или иные методы измерения.
        Если умножить погрешность на 1/10 погрешность уменьшиться, а если на 0 так вообще пропадёт. Но какой толк от этих операций? Нам же нужна сама величина и доверительный интервал. Какой смысл «заверать» ошибку измерения?
          0
          Как известно, погрешность имеет две составляющие: систематическую и случайную, фильтрация уменьшает случайную составляющую, именно для этого и применяется фильтр Калмана. Энтропийная форма представления погрешности отображает интервал неопределённости (не путайте с доверительным интервалом). По Вашему мнению, применение фильтра не снижает случайную составляющую погрешности, обусловленную наличием шумов — Вы это серьёзно?!.. Разберитесь с терминологией…
            –1
            Каких шумов? У вас есть шум в виде наводок 50Гц? У вас есть физическая модель системы — константа. Т.е. грубо вы линейкой измеряете длину 55 раз. Систематическая погрешность равна разрешающей способности линейки, а случайная зависит от факторов которые вам не подвластным (влажность, температура, внутренняя структура и т.п.). Вы оцениваете дисперсию полученных значений, умножаете на коэффициент (определяющий вероятность попасть в доверительный интервал) и на основе этих цифр получаете доверительный интервал.

            «Энтропийная форма представления погрешности отображает интервал неопределённости» — откуда и куда отображает?

            Фильтры есть разные и они предназначены для подавления разного рода помех. Но для этого надо иметь модель помех или модель процесса.
              +1
              Фильтр Калмана в каждом цикле корректирует дисперсию, и, при многократных измерениях, дисперсия случайной величины уменьшается, а следовательно сужается интервал неопределённости и уменьшается энтропийная погрешность. Закон распределения случайной погрешности определяется по энтропийному коэффициенту и контрэксцессу и не зависит от физической природы случайной величины.
                0
                Вы не можете менять дисперсию случайной величины. Вы можете только оценить её величину (закон распределения).
                  +1
                  «Вы не можете менять дисперсию случайной величины. Вы можете только оценить её величину (закон распределения).»
                  Да нужно говорить об оценке дисперсии для ограниченной выборки по которой мы её определяем, но это не как не влияет на факт уменьшения оценки дисперсии при использовании фильтра Калмана в режиме многократных измерений.
                    0
                    Усреднение результатов многократных наблюдений при постоянстве значения измеряемой величины является наиболее эффективным методом уменьшения случайной погрешности измерения. При проведении многократных (n) наблюдений одного и того же значения физической величины во многих случаях в качестве результата измерения выбирается среднее значение результатов наблюдений. В этом случае среднее квадратическое отклонение результата измерения уменьшается в n раз.
                      0
                      Уточняю: В этом случае среднее квадратическое отклонение результата измерения уменьшается в корень квадратный из n
                        0
                        А вы проверьте:
                        import scipy.stats
                        from random import normalvariate
                        
                        class Disp:
                            def __init__(self):
                                self.m=0
                                self.d=0
                                self.n=0
                            def add(self,x):
                                self.n=self.n+1
                                if self.n>1:
                                    self.d=self.d*(self.n-2)/(self.n-1)+(x-self.m)**2/self.n
                                self.m=self.m+(x-self.m)/self.n
                            def getN(self):
                                return self.n # =n
                            def getM(self):
                                return self.m # =sum(x[i])/n
                            def getD(self): 
                                return self.d # =sum((x[i]-m)^2)/(n-1)
                        
                        def st(alpha,n):
                            return scipy.stats.t.ppf((1+alpha)/2, n-1)
                        
                        disp=Disp()
                        xo=1
                        mu=0.5
                        print("Estimate xsi(xo=%.2f,mu=%.2f)\n"%(xo,mu))
                        alpha=0.95
                        print("  count     xsi    mean    disp  err(%.2f)"%alpha)
                        j=10
                        for i in range(100000):
                            y=normalvariate(xo,mu)
                            disp.add(y)
                            if i+1==j:
                                j=j*10
                                n,m,d=disp.getN(),disp.getM(),disp.getD()**0.5
                                print("%7d\t%7.4f\t%7.4f\t%7.4f\t%7.4f"%( n, y, m, d, d*st(alpha,n) ))
                        

                        Estimate xsi(xo=1.00,mu=0.50)
                        
                          count     xsi    mean    disp  err(0.95)
                             10	 0.4935	 0.9605	 0.5443	 1.2313
                            100	 0.8873	 1.0034	 0.4685	 0.9296
                           1000	 1.3086	 0.9905	 0.4961	 0.9734
                          10000	 1.0788	 1.0022	 0.5023	 0.9845
                         100000	 0.8025	 0.9999	 0.5000	 0.9800

                        Оценка дисперсии не зависит от n, она сходится к значению дисперсии случайной величины, но никак не уменьшается.
        0
        2019 год на дворе, а всё старика Калмана лохматят.
          0
          Проверяю!
          from numpy import *
          import matplotlib.pyplot as plt
          from scipy.stats import norm
          n_iter = 100 # Число итераций.
          sz = (n_iter,) # Размер массива
          x =2# Истинное значение измеряемой величины (фильтру неизвестно)
          R1 = 0.1 # Ср. кв. ошибка измерения.
          R = R1*R1 # Дисперсия
          nr=«нормальным распределением»
          y=norm.rvs( x, R1, size=sz)
          Q = 1e-5 # Дисперсия случайной величины в модели системы
          # Выделение памяти под массивы:
          xest1 = zeros(sz) # Априорная оценка состояния
          xest2 = zeros(sz) # Апостериорная оценка состояния
          P1 = zeros(sz) # Априорная оценка ошибки
          P2 = zeros(sz) # Апостериорная оценка ошибки
          G = zeros(sz) # Коэффициент усиления фильтра
          xest2[0] = 0.0
          P2[0] = 1.0
          for k in arange(1, n_iter,1): # Цикл по отсчётам времени.
          xest1[k] = xest2[k-1] # Априорная оценка состояния.
          P1[k] = P2[k-1] + Q# Априорная оценка ошибки.
          # После получения нового значения измерения вычисляем апостериорные оценки:
          G[k] = P1[k] / ( P1[k] + R )
          xest2[k] = xest1[k] + G[k] * ( y[k] — xest1[k] )
          P2[k] = (1 — G[k]) * P1[k]
          plt.title('Ошибки при подавлении шумов \n с %s'%nr, size=12)
          valid_iter = arange(1, n_iter,1) # P1 на 0 м шаге не определено
          plt.plot(valid_iter, P1[valid_iter])
          plt.xlabel('Номер итерации')
          plt.ylabel('Априорная оценка ошибки')
          plt.setp(plt.gca(), 'ylim', [0, .01] )
          plt.show()
            0
            А если отложить модуль ошибки пунктиром на каждой итерации и сравнить с оценкой ошибки?
            image
            Вас не смущает оценка ошибки и сама ошибка?
            code.py
            from numpy import *
            import matplotlib.pyplot as plt
            from scipy.stats import norm
            
            n_iter = 1000 # Число итераций.
            sz = (n_iter,) # Размер массива
            x =2 # Истинное значение измеряемой величины (фильтру неизвестно)
            R1 = 0.1 # Ср. кв. ошибка измерения.
            R = R1*R1 # Дисперсия
            nr="нормальным распределением"
            y=norm.rvs( x, R1, size=sz)
            Q = 1e-5 # Дисперсия случайной величины в модели системы
            xest1 = zeros(sz) # Априорная оценка состояния
            xest2 = zeros(sz) # Апостериорная оценка состояния
            P1 = zeros(sz)    # Априорная оценка ошибки
            P2 = zeros(sz)    # Апостериорная оценка ошибки
            G = zeros(sz)     # Коэффициент усиления фильтра
            err1=zeros(sz)    # Ошибка оценки состояния
            xest2[0] = 0.0
            P2[0] = 1.0
            for k in arange(1, n_iter,1): # Цикл по отсчётам времени.
                xest1[k] = xest2[k-1]     # Априорная оценка состояния.
                P1[k] = P2[k-1] + Q       # Априорная оценка ошибки.
                # После получения нового значения измерения вычисляем апостериорные оценки:
                G[k] = P1[k] / ( P1[k] + R )
                xest2[k] = xest1[k] + G[k] * ( y[k] - xest1[k] )
                P2[k] = (1 - G[k]) * P1[k]
                err1[k]=abs(xest1[k]-x)
            
            plt.title('Ошибки при подавлении шумов \n с %s'%nr, size=12)
            i = arange(1, n_iter,1) # P1 на 0 м шаге не определено
            plt.plot(i, P1[i])
            plt.plot(i, err1[i],'--')
            plt.xlabel('Номер итерации')
            plt.ylabel('Априорная оценка ошибки')
            plt.setp(plt.gca(), 'ylim', [0, 0.05] )
            plt.show()
            

              0
              Фильтр имеет свой алгоритм работы. То что Вы сделали отношения к ошибке не имеет а свидетельствует только о неизменности дисперсии шума !!!..
                0
                О терминологии- Оценка погрешности численного интегрирования. Различают два вида оценок априорные и апостериорные. Априорную оценку получают заранее, до проведения расчетов, на основе теоретического анализа квадратурной формулы. Апостериорную оценку определяют после вычислений на основе сопоставления результатов расчетов, проведенных при разных числах отрезков разбиения
                  –1
                  А причем тут численное интегрирование. Априорную оценку мы знаем т.к. сами задали распределение случайной величины. А апостериорную можно определить оценив результат с полученный с разной степенью точности. Для численного интегрирования есть простой метод оценки — процесс Эйткена.
                  Но какое отношение это имеет к вашим измерениям?

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

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

              Самое читаемое