18402218@smail.cczu.edu.com
感谢博文:
https://blog.csdn.net/zhangyue_lala/article/details/64437547?depth_1-utm_source=distribute.pc_relevant.none-task&utm_source=distribute.pc_relevant.none-task
基于牛顿插值的实时概率预测方法
前言:实时预测作用于采样频率较高的场景(大于10Hz),实时预测的作用大多是为了提高采样频率,类似于提前插值。从实际上说,未来时刻的位置信息是未知的,只能基于概率预测。
注意!!!预测前的数据必须平滑化处理
1.实时预测的原则
1.尽可能多地利用确定的坐标信息(Newton法)
2.更多地相信时间近的坐标信息(定最小确性度)
3.具有安全性(长时间没有得到实际坐标点的更新时的安全处理,概率监测)
2.算法
2.1对于n个已知点的概率预测
已知n个x1,x2,x3…xn
下标数字越大时间离得越近
只有1个点xn时:
当
Δ
\Delta
ΔT–>0时,x在xn位置的概率趋近于1,即:
P(xn)–>1,其它地方的概率均为0。
在
Δ
\Delta
ΔT时间上我们可以知道x的位置概率服从高斯分布:
N
N
N(xn,
p
2
p^2
p2)
当
Δ
\Delta
ΔT–>0时,p–>
0
0
0,为了计算p取0.1
随着时间的增长,x位置的不确定性也随之增加,对于高斯分布,即p值会慢慢变大,直到平坦。(图形如下)
这里
z
z
z=
N
N
N(0,
p
2
p^2
p2)
p
p
p=y
这里我们让方差
p
p
p随着时间线性变化,随着时间的变化,虽然0点的位置概率始终最大,但概率值在不断缩小,当
p
p
p增长到一定程度,0点的可行度是很低的,我们不得不放弃这个点,所以在模型中我们必须给出最小确信度
P
P
Pmin用来舍弃一些时间较长的点。
如果最相近的一点低于最小确信度,进入报错程序
只有2个点x(n-1),xn时:
当
Δ
\Delta
ΔT–>0时,x在xn位置的概率趋近于1,即:
P(xn)–>1,其它地方的概率均为0。
如果x(n-1)点在我们上述的最小确性度内:
我们可以计算出
V
V
V=[
f
f
f(xn)-
f
f
f(x(n-1))]/xn-x(n-1)
在
Δ
\Delta
ΔT时间上我们可以知道x的位置概率服从高斯分布:
N
N
N(xn+
V
V
V*
Δ
\Delta
ΔT,
p
2
p^2
p2)
这时图像如下:
这里
z
z
z=
N
N
N(
p
/
2
p/2
p/2,
p
2
p^2
p2)
p
p
p=y
这里我们让方差
p
p
p随着时间线性变化,随着时间的变化,虽然p/2点的位置概率始终最大,但概率值在不断缩小,当
p
p
p增长到一定程度,当x(n-1)点已经小于我们的最小确性度
P
P
Pmin,我们就要使用上面的只有一个点的预测方法。
同理有n个已知点可以用上述方法预测
以上方法满足了我们的2,3条原则:
2.更多地相信时间近的坐标信息(定义最小确性度)
3.具有安全性(长时间没有得到实际坐标点的更新时的安全处理,概率监测)
2.2定义差商
可以很好地得到这些性质:
定义差商对我们的下面表示方便了不少。
可以看到:
对于物体运动,一阶差商得到速度,二阶差商得到加速度,三阶差商是加速度的微分…对于n阶差商我们会用到n+1个状态的坐标。对于下一个
Δ
\Delta
ΔT时刻的位置信息的预测,我们可以从n阶差商逐渐下降到0阶差商得到下一个
Δ
\Delta
ΔT时刻的坐标:
Δ
P
t
\Delta Pt
ΔPt
=
f
f
f[x(n-1,xn)]+
f
f
f(xn)
={
f
f
f[x(n-2),x(n-1),xn]
Δ
\Delta
ΔT+
f
f
f[x(n-1,xn)]}
Δ
\Delta
ΔT+
f
f
f(xn)
…
通过牛顿法,我们可以达到我们的第一条原则:
1.尽可能多地利用确定的坐标信息
2.3算法流程
3.python实现
3.1加上噪音的模拟数据
import math
import json
import numpy
import threading
import time
import matplotlib.pylab as plt
#store
pool=[]
#store predicted position
predict_1=[]
#store predicted time
time_1=[]
#threading protected lock
lock=threading.Lock()
#to get Difference quotient of a list
def get_DQ(list1=[]):
answers=[]
next=list1
new_next=[]
while len(next)>1:
new_next=[]
for i in range(len(next)-1):
temp=[]
val=(next[i+1][0]-next[i][0])/(next[i+1][1]-next[i][1])
temp.append(val)
time1=(next[i+1][1]+next[i][1])/2
temp.append(time1)
new_next.append(temp)
if i==len(next)-2:
answers.append(round(val,4))
next=new_next
return answers
#class a class of predict ,with inherited class of threading,Thread
class predict(threading.Thread):
#unit of od trust:microseconds,of freq:Hz,
def __init__(self,trust,dimension,freq):
threading.Thread.__init__(self)
self.trust=trust
self.dimension=dimension
self.freq=freq
self.piece=1000/freq
self.thisTime=0
self.last_time=0.0
self.DQ=[]
self.next_predict=0
self.last_position=0
def run(self):
global pool
global predict_1
global lock
while True:
lock.acquire()
for i in pool:
if i[1]<self.thisTime-self.trust:
pool.pop(0)
else :
break
lock.release()
length=len(pool)
if len(pool)==0:
self.go_error()
elif pool[length-1][1]!=self.last_time:
predict_1.append(pool[length-1][0])
self.last_time=pool[length-1][1]
self.last_position=pool[length-1][0]
self.thisTime=self.last_time
time_1.append(self.thisTime)
else :
self.DQ=get_DQ(pool)
value=self.calculate()
self.last_position+=value
predict_1.append(self.last_position)
time_1.append(self.thisTime)
time.sleep(self.piece/1000)
self.thisTime+=self.piece
def go_error(self):
print("out of time!")
if len(predict_1)>10:
exit(0)
def calculate(self):
predict_2=0
size=len(self.DQ)
for i in range(len(self.DQ)):
predict_2=(predict_2+self.DQ[size-1-i])*self.piece
return predict_2
def main():
#list1=[[1,2.3],[2.1,3],[2.9,4],[4,5.3],[5.1,6],[6.2,7.3],[7.2,9],[8,10],[8.8,11],[9.7,12]]
#list2=[[1,2],[2,3],[3,4],[4,5],[5,6],[6,7],[7,8],[7.3,9],[7.4,10]]
#to use this class with trust =200,dimension=1,expect freq=100
task1=predict(200,1,100)
global pool
global lock
#read json file
with open("myData3.json",'r') as load_f:
load_dict = json.load(load_f)
#threading start
task1.start()
last_time=0
#store original data
or_p=[]
or_t=[]
for i in load_dict:
time.sleep((load_dict[str(i)][1]-last_time)/1000)
lock.acquire()
pool.append(load_dict[str(i)])
lock.release()
or_p.append(load_dict[str(i)][0])
or_t.append(load_dict[str(i)][1])
last_time=load_dict[str(i)][1]
task1.join()
#draw the graph
plt.plot(time_1,predict_1,'.')
plt.plot(or_t,or_p,'.')
plt.plot()
plt.show()
if __name__=="__main__":
main()
实验中用到的数据可以从这里下载
https://www.jianguoyun.com/p/DVC73JYQ7L2UCBiHh-4C
代码实时地预测,从25帧达到100帧,图像如下:
图像中蓝色点为插值插入的点,橙色为采样数据点。
我们可以得出:
1.该预测可以很好地预测出足够平滑采样点的数据
(要求采样帧率至少10Hz)
2.该预测是非线性的,从最后的一段蓝色的预测中可以得出
3.如果采样数据不够平滑,该预测结果离散度将会很大
无噪音的函数图像数据
我们再对正弦函数做一次预测,以计算一次偏差率
可以得到这样的模拟图像:
在这次模拟中,采样频率为12.5Hz,期待帧率为100帧,可以发现出现了明显的滞后现象。
采样频率设置为25Hz后,滞后现象有明显改进:
最后我们来求出12.5Hz时的累积误差,公式如下:
D
D
Dsum=
∑
\sum
∑abs(
f
f
f(x)-
P
r
e
d
i
c
t
Predict
Predict(x))
最后累积误差小于4.0
本人不擅数理,欢迎批评指正