Hallo Filter

1
2
3
4
5
6
7
import numpy as np
from scipy import stats
import matplotlib.pyplot as plt
# %matplotlib ipympl
%matplotlib widget

figsize=(20, 5)

噪声的数学模型

先来回想一下梦开始的地方，小学数学里测量一个物体的长度、重量巴拉巴拉等物理量时，为了减小误差的方法是————多次测量取平均值。
所以多次测量取平均值是怎么减小误差的呢？或者用数学语言来描述的话，对于未知随机变量 XX ，进行多次实验后，其均值 Xˉ\bar X 的抽样分布如何变化？

噪声的数学描述

中心极限定理告诉我们: 只要随机变量相互独立，每个随机变量对和的影响都是微小的，哪怕它们的分布类型不同，其和标准化后都有标准正态的极限分布。话句话说，某一事件受到许多相互独立且微小的随机因素影响，总的影响就可以看作服从或近似服从正态分布。
不妨来仿真一下，下述代码测试计算了 50000 次随机事件，每次随机事件均为 100 次平均分布的随机变量之和。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
import math

# 绘制标准正态分布曲线
μ = 0
σ = math.sqrt(1)
x = np.linspace(μ - 3*σ, μ + 3*σ)
y = np.exp(-(x - μ) ** 2 / (2 * σ ** 2)) / (math.sqrt(2*math.pi)*σ)
plt.figure(figsize=figsize)
plt.plot(x, y)

# 中心极限定理仿真
a = []
for i in range(50000):
    a.append(np.sum(np.random.uniform(10, 20, 100)))
a = (a - np.average(a)) / np.std(a)
plt.hist(a, 100, density=True)
Figure

非常amazing啊！100次随机试验 U(10,20)U(10,20) 经过标准化后的的频率分布直方图竟然被标准正态分布函数完美的包络了。这就是正态分布在生活中如此常见的原因。
继续以测量误差为例，其受到许多相互独立随机因素的影响，如测量环境温度、湿度，测量工具的精密程度以及测量者的心理因素、测量的态度等的影响，而每种影响都不占主要地位，那么它们总和造成的总误差就近似地服从正态分布。

单正态总体的抽样分布

现在我们知道误差在数学上的模型了，但多次测量为什么能减小误差呢？
由上一小点易知，对于某一可直接观测物理量的单次试验 E，设其测量结果为随机事件 XX ，且 X∼N(μ,σ2)X\sim N(μ,σ^2)。
而我们的目的是求出未知参量 μμ ,而为了求出未知参量，需要观测。设进行了 n 次相互独立的观测，其结果分别为 XiX_i 。

Xˉ=1n∑i=1nXi∼N(μ,1nσ2)\bar X=\frac{1}{n}\sum_{i=1}^{n}X_i \sim N(\mu, \frac{1}{n}\sigma^2)

不难看出，方差与次数成反比。所以，随测量次数的增加方差会快速下降。如果你觉得不够直观的话不妨用概率来描述，由切比雪夫不等式

P{∣Xˉ−E(Xˉ)∣≥ε}≤D(Xˉ)ε2==>P{∣Xˉ−μ∣≥ε}≤σ2nε2P\{|\bar X - E(\bar X)| \geq \varepsilon\} \leq \frac{D(\bar X)}{\varepsilon^{2}} \qquad ==> \qquad P\{|\bar X - \mu | \geq \varepsilon\} \leq \frac{\sigma ^2}{n\varepsilon^{2}}

可知，事件 “随机变量 Xˉ\bar X 关于其期望产生了偏差 ε” 发生的概率上限为 σ2nε2\frac{\sigma ^2}{n\varepsilon^{2}} 。不妨令 n=k/ε2n = k/ε^2 于是有

P{∣Xˉ−E(Xˉ)∣≥ε}≤σ2kP\{|\bar X - E(\bar X)| \geq \varepsilon\} \leq \frac{\sigma^2}{k}

这就是说，当实验次数为 kε2\frac{k}{ε^2} 时，随机变量 Xˉ\bar X 关于其期望产生了偏差 ε 的概率上限为 σ2k\frac{\sigma^2}{k}

从直观感觉上来说的话，每次测量都包含了“信息”。每一次新的测量都代表了“新息”的进入，多次测量使得系统的信息增多，不确定性下降，从而使得新息的信息熵降低。在这个角度上倒是很容易理解为什么多次测量的边际效应会十分显著了。

滤波器初体验

滤波器仿真

对于一个静态系统来说，也没什么花里胡哨的玩法，暴力求平均即可。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
### 模型参数 ###
# 采样参数
SampleFre = 10 
TotalTime = 10 
Step = 1/SampleFre 
SampleNum = SampleFre * TotalTime 
#观测序列
RealValue = []        #实际值
MeasValue = []       #测量值
EstiSequence = []        #滤波值

### 仿真 ###
x = 0
while x < TotalTime:
    μ = 1                                #实际值
    RealValue.append(μ)

    temp = np.random.normal(μ, 1)      #测量值 
    MeasValue.append(temp)

    EstiSequence.append(np.average(MeasValue)) #滤波值

    x += Step

### 绘图展示 ###
plt.figure(figsize=figsize)
plt.plot(RealValue, 	'*-', label="Actual", markersize=5)
plt.plot(MeasValue, label="Measure", markersize=5)
plt.plot(EstiSequence, 	'.-', label="Filter", markersize=5)
plt.legend(loc='lower right') #绘制图例
Figure

emmmmmm 效果不错…吧？
随着次数的增多"感觉"误差越来越小了，但该怎么用数学语言描述呢？

  1. 事物的好坏、多少等等，都是定性分析，是有潜在参考标准的，单独讨论没有任何意义。作为工程师，我们需要将定性分析转化为定量分析，具体的去分析某种事物的性能。
    记得高中学集合时老师举的一个例子：“把所有高个子的人记作一个集合”，这种说法是错误的。因为高是一个相对词，正确的说法应该是所有身高大于 180cm 的人记作一个集合。
  2. 任何改进都一定会引入新的缺点，但往往我们只关注某一方面的性能，在另一方面上的缺点是可以被忽略的。所以在工程上千万不要有强者恒强的想法，而是要根据情况具体分析后选择最优的解法。

滤波器性能的分析

现在我们通过 n 次观测得到了 n 个样本的均值，记作

x^n=xˉn\hat{x}_{n}=\bar{x}_{n}

即对 μ\mu 的点估计。

那么这个估计有多可靠呢？想要判断估计值的可靠性，我们需要找到估计值 x^n\hat{x}_{n} 与真值 μ\mu 之间的关系，用数学语言表达：试求一函数 G=G(x^n,μ)G=G(\hat{x}_{n}, \mu)， 并求出 GG 的分布密度。

即找到一个枢轴量，并求该枢轴量的分布。一般的对于正态总体，其枢轴量为常见分布，如正态分布、卡方分布、t 分布、F 分布。

让我们再次回到最开始的地方，

Xˉ=1n∑i=1nXi∼N(μ,1nσ2)\bar X=\frac{1}{n}\sum_{i=1}^{n}X_i \sim N(\mu, \frac{1}{n}\sigma^2)

易得到

Xˉ=ΔX^∼N(μ,1nσ2)\bar X\overset{\Delta}{=} \hat X \sim N(\mu, \frac{1}{n}\sigma^2)

你就不想对它标准化一下么？

G=X^−μσ/nG = \frac{\hat X - \mu} { \sigma / \sqrt{n} }

测量方差已知的情况

当方差 σ\sigma 已知的时候，枢轴量 GG 的分布显然满足 G N(0,1)G~N(0,1) 。
此时， 95% 对称双侧置信区间为

(Xˉ−σnZ2.5%,Xˉ+σnZ2.5%)\left(\bar{X}-\frac{\sigma}{\sqrt{n}} Z _{2.5\%}, \bar{X}+\frac{\sigma}{\sqrt{n}} Z_{2.5\%} \right)

Z分布即标准正态分布。毕竟不是数学教程，这里就不详细展开置信区间了。学过数一都应该会的。
简单来说，95% 置信区间的意思是“我们有 95% 的信心说这个置信区间包含了真值。即如果大量重复试验并计算置信区间，那么应该有 95% 的置信区间包含这个真值”。这与某个置信区间有 95% 的概率包含了估计的真值是不一样的。

测量方差未知的情况

当方差 σ\sigma 未知的时候，使用样本标准差 SS 代替总体标准差 σ\sigma，此时枢轴量 GG 的分布满足 Gt˜(n−1)G\~t(n-1) 。
此时， 95% 对称双侧置信区间为

(Xˉ−σntα/2(n−1),Xˉ+σntα/2(n−1))\left(\bar{X}-\frac{\sigma}{\sqrt{n}} t_{\alpha / 2}(n-1), \bar{X}+\frac{\sigma}{\sqrt{n}} t_{\alpha / 2}(n-1)\right)

当 n 充分大的时候(n>30n>30)，一般也可以直接使用标准正态分布来计算。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
#观测序列
RealSequence = []       #实际值
err0 = []               #误差
std0 = []               #标准差
MeasSequence = []       #测量值
err1 = []               #误差
std1 = []               #标准差
EstiSequence = []       #滤波值

### 仿真 ###
x = 0
while x < TotalTime:
    trueValue = 1
    RealSequence.append(trueValue)

    temp = np.random.normal(trueValue, 10)		#测量值 
    MeasSequence.append(temp)
    err0.append(np.abs(temp - trueValue))		#误差

    temp = np.average(MeasSequence)				#滤波值
    EstiSequence.append(temp)
    std0.append(np.std(MeasSequence)/np.sqrt(np.size(MeasSequence)) * stats.t.isf(0.025, np.size(err0)-1))			# 95%置信区间

    err1.append(np.abs(temp - RealSequence[-1]))   #误差
    std1.append(np.std(MeasSequence)/np.sqrt(np.size(MeasSequence)) * stats.t.isf(0.025, np.size(err1)-1))			# 95%置信区间

    x += Step
1
2
3
4
5
6
7
8
9
10
11
12
### 绘图展示 ###
plt.figure(figsize=figsize)
plt.plot(RealSequence, 	'*-', label="Actual", markersize=5)
plt.plot(MeasSequence, 	'o-', label="Measure", markersize=5)
plt.errorbar(range(0, np.size(EstiSequence)), EstiSequence, yerr=std0, label="Filter", fmt='.-', capthick=2, capsize=10)
plt.legend(loc='lower right') #绘制图例

plt.figure(figsize=figsize)
plt.plot(err0, '*-', label="Measure Error")
plt.plot(err1, '.-', label="Estimate Error")
plt.plot([0 for i in range(0, np.size(err1))], '-', markersize=5)
plt.legend(loc='lower right') #绘制图例
Figure
Figure

多运行几次后，大致可以总结成：

  • 前期(n≈30n\approx 30)一般收敛较快。
  • 误差的绝对值并不总是一直减小的。

均值的递推公式

xˉ=∑i=1nxinxˉ=∑i=1n−1xin−1⋅n−1n+xnnxˉ=xˉlast⋅(n−1)+xnnxˉ=xˉlast+1n(xn−xˉlast)\begin{align} \nonumber \bar{x} & = \frac{\sum_{i=1}^{n}{x_i}}{n}\\ \nonumber \bar{x} & =\frac{\sum_{i=1}^{n-1}{x_i}}{n-1} \cdot \frac{n-1}{n} + \frac{x_n}{n}\\ \nonumber \bar{x} & = \frac{\bar{x}_{last} \cdot (n-1) + x_n}{n}\\ \nonumber \bar{x} & = \bar{x}_{last} + \frac{1}{n}(x_n-\bar{x}_{last}) \end{align}

上述形式的在应用中的优点主要是:

  • 不需要储存全部的测量值，只需要保存当前次数与平均值即可(int + float)。
  • 计算次数少，最后的形式只需计算两次加法与一个除法。

当然，前文提到过任何改进都将引入新的缺点：

  • 计算机的float类型是不可靠的，每次迭代都会带来误差。
  • 与最开始小学生都会的形式相比递推形式更难了一些。(啥？你说这点难度也算缺点？记住这个形式，以后你还会见到他的🤭)