#下面的代码也是完全基于上述作者的库函数完成的，所以需要先去Github下载库函数安装,或者直接使用
pip install filterpy
#但是作者忘记把离散的白噪声函数放进库里，那段函数可以在书里找到，或者改天我可以放在网盘里分享出来
#不影响下面代码展示

#离散贝叶斯滤波器的实现展示
#The Kalman filter belongs to a family of filters called Bayesian filters.
#NMEFC YuHao
import numpy as np
import matplotlib.pyplot as plt
from filterpy.discrete_bayes import normalize
from filterpy.discrete_bayes import predict
from filterpy.discrete_bayes import updatex = range(10)
hallway = np.array([1,1,0,0,0,0,0,0,1,0])
#Tracking a Dogdef update_belief(hall,belief,z,correct_scale):for i ,val in enumerate(hall):if val == z:belief[i] *= correct_scale
belief = np.array([0.1]*10)
reading = 1
update_belief(hallway,belief,z = reading,correct_scale = 3.)
print("belief:",belief)
print("sum = ",sum(belief))
plt.figure()
plt.bar(x,belief)
plt.title("belief")
plt.show()#归一化
homogeneity_array = belief/sum(belief)
print(homogeneity_array)
print(hallway == 1 )
#计算后验概率
def scale_update(hall ,belief, z, z_prob):scale = z_prob / (1. - z_prob)belief[hall == z] *= scalenormalize(belief)
belief = np.array([0.1]*10)
scale_update(hallway ,belief ,z=1,z_prob = .75)
print('sum =',sum(belief))
print("probability of door =",belief[0])
print("probability of wall =", belief[2])
plt.figure()
plt.bar(x,belief)
plt.title("belief")
plt.show()def scale_update(hall,belief,z,z_prob):scale = z_prob / (1. - z_prob)likelihood = np.ones(len(hall))likelihood[hall == z] *= scalereturn normalize(likelihood * belief)
def lh_hallway(hall, z, z_prob):#compute likelihood that a measurement matches#positions in the hallway.try:scale = z_prob / (1. - z_prob)except ZeroDivisionError:scale = 1e8likelihood = np.ones(len(hall))likelihood[hall==z] *= scalereturn likelihoodbelief = np.array([0.1] * 10)
likelihood = lh_hallway(hallway, z=1, z_prob=.75)
print(normalize(likelihood * belief))
def perfect_predict(belief,move):n = len(belief)result = np.zeros(n)for i in range(n):result[i] = belief[(i-move) % n]return result
belief = np.array([.35,.1,.2,.3,0,0,0,0,0,.05])belief = perfect_predict(belief,1)
plt.bar(x,belief)
plt.title("belief")
plt.show()
#添加噪声def predict_move(belief, move, p_under, p_correct, p_over):n = len(belief)prior = np.zeros(n)for i in range(n):prior[i] = (belief[(i-move) % n]   * p_correct +belief[(i-move-1) % n] * p_over +belief[(i-move+1) % n] * p_under)      return priorbelief = [0., 0., 0., 1., 0., 0., 0., 0., 0., 0.]
prior = predict_move(belief, 2, .1, .8, .1)
print(prior)
fig = plt.figure()
ax1 = fig.add_subplot(1,2,1)
ax2 = fig.add_subplot(1,2,2)
ax1.bar(x,belief)
ax2.bar(x,prior)
ax1.set_title("belief")
ax2.set_title("prior")
ax1.set_xticks(x)
ax2.set_xticks(x)
plt.show()
belief = [0, 0, .4, .6, 0, 0, 0, 0, 0, 0]
prior = predict_move(belief, 2, .1, .8, .1)
fig = plt.figure()
ax1 = fig.add_subplot(1,2,1)
ax2 = fig.add_subplot(1,2,2)
ax1.bar(x,belief)
ax2.bar(x,prior)
ax1.set_title("belief")
ax2.set_title("prior")
ax1.set_xticks(x)
ax2.set_xticks(x)
plt.show()belief = np.array([1.0,0,0,0,0,0,0,0,0,0])
beliefs = []
for i in range(100):belief = predict_move(belief,1,.1,.8,.1)beliefs.append(belief)
print("Final Belief:",belief)
print(len(beliefs))
print(beliefs[1])
for i in range(len(beliefs)):plt.cla()plt.bar(x,beliefs[i])plt.xticks(x)plt.pause(0.05)plt.show()
"""
#卷积泛化
#卷积算法的python描述，但是时间复杂度太高，pass
def predict_move_convolution(pdf ,offset,kernel):N = len(pdf)kN = len(kernel)width = int((kN - 1)/2)prior = np.zeros(N)for i in range(N):for k in range(kN):index = (i + (width - k)-offset)%Nprior[i] += pdf[index] * kernel[k]return prior
"""belief = [.05, .05, .05, .05, .55, .05, .05, .05, .05, .05]
prior = predict(belief,offset = 1,kernel = [.1,.8,.1])
fig = plt.figure()
ax1 = fig.add_subplot(1,2,1)
ax2 = fig.add_subplot(1,2,2)
ax1.bar(x,belief)
ax2.bar(x,prior)
ax1.set_title("belief")
ax2.set_title("prior")
ax1.set_xticks(x)
ax2.set_xticks(x)
plt.show()
prior = predict(belief, offset=3, kernel=[.05, .05, .6, .2, .1])
#make sure that it shifts the positions correctly for movements greater than one and for asymmetric kernels
#making a prediction of where the dog is moving, and convolving the probabilities to get the priorfig = plt.figure()
ax1 = fig.add_subplot(1,2,1)
ax2 = fig.add_subplot(1,2,2)
ax1.bar(x,belief)
ax2.bar(x,prior)
ax1.set_title("belief")
ax2.set_title("prior")
ax1.set_xticks(x)
ax2.set_xticks(x)
plt.show()
"""
#Integrating Measurements and Movement Updates
"""
"""
请仔细阅读下一段话
Let's think about this intuitively.Consider a simple case - you are tracking a dog while he sits still.During each prediction you predict he doesn't move.Your filter quickly converges on an accurate estimate of his position.Then the microwave in the kitchen turns on, and he goes streaking off.You don't know this, so at the next prediction you predict he is in the same spot.But the measurements tell a different story.As you incorporate the measurements your belief will be smeared along the hallway, leading towards the kitchen.On every epoch (cycle) your belief that he is sitting still will get smaller,and your belief that he is inbound towards the kitchen at a startling rate of speed increases.
最精彩的部分来了：反馈器"""
tank1 = []
tank2 = []hallway = np.array([1, 1, 0, 0, 0, 0, 0, 0, 1, 0])
prior = np.array([.1] * 10)
likelihood = lh_hallway(hallway, z=1, z_prob=.75)
posterior = normalize(likelihood * prior)
tank2.append(posterior)
tank1.append(prior)
fig = plt.figure()
ax1 = fig.add_subplot(1,2,1)
ax2 = fig.add_subplot(1,2,2)
ax1.bar(x,prior)
ax2.bar(x,posterior)
ax1.set_title("prior")
ax2.set_title("posterior")
ax1.set_xticks(x)
ax1.set_ylim(0,0.5)
ax2.set_ylim(0,0.5)
ax2.set_xticks(x)
plt.show()#After the first update we have assigned a high probability to each door position, and a low probability to each wall position.
prior = np.array([.1] * 10)
kernel = (.1, .8, .1)
prior = predict(posterior, 1, kernel)
tank1.append(posterior)
tank2.append(prior)
fig = plt.figure()
ax1 = fig.add_subplot(1,2,1)
ax2 = fig.add_subplot(1,2,2)
ax1.bar(x,prior)
ax2.bar(x,posterior)
ax1.set_title("posterior")
ax2.set_title("prior")
ax1.set_xticks(x)
ax1.set_ylim(0,0.5)
ax2.set_ylim(0,0.5)
ax2.set_xticks(x)
plt.show()
#The predict step shifted these probabilities to the right, smearing them about a bit. Now let's look at what happens at the next sense.likelihood = lh_hallway(hallway, z=1, z_prob=.75)
posterior = update(likelihood, prior)
tank1.append(prior)
tank2.append(posterior)
fig = plt.figure()
ax1 = fig.add_subplot(1,2,1)
ax2 = fig.add_subplot(1,2,2)
ax1.bar(x,prior)
ax2.bar(x,posterior)
ax1.set_title("prior")
ax2.set_title("posterior")
ax1.set_xticks(x)
ax1.set_ylim(0,0.5)
ax2.set_ylim(0,0.5)
ax2.set_xticks(x)
plt.show()
#Notice the tall bar at position 1. This corresponds with the (correct) case of starting at position 0, sensing a door, shifting 1 to the right, and sensing another door. No other positions make this set of observations as likely. Now we will add an update and then sense the wall.
prior = predict(posterior, 1, kernel)
likelihood = lh_hallway(hallway, z=0, z_prob=.75)
posterior = update(likelihood, prior)
tank1.append(prior)
tank2.append(posterior)
fig = plt.figure()
ax1 = fig.add_subplot(1,2,1)
ax2 = fig.add_subplot(1,2,2)
ax1.bar(x,prior)
ax2.bar(x,posterior)
ax1.set_title("prior")
ax2.set_title("posterior")
ax1.set_xticks(x)
ax1.set_ylim(0,0.5)
ax2.set_ylim(0,0.5)
ax2.set_xticks(x)
plt.show()
#This is exciting! We have a very prominent bar at position 2 with a value of around 35%. It is over twice the value of any other bar in the plot, and is about 4% larger than our last plot, where the tallest bar was around 31%. Let's see one more cycle.
prior = predict(posterior, 1, kernel)
likelihood = lh_hallway(hallway, z=0, z_prob=.75)
posterior = update(likelihood, prior)
prior = predict(posterior, 1, kernel)
likelihood = lh_hallway(hallway, z=0, z_prob=.75)
posterior = update(likelihood, prior)
tank1.append(prior)
tank2.append(posterior)
fig = plt.figure()
ax1 = fig.add_subplot(1,2,1)
ax2 = fig.add_subplot(1,2,2)
ax1.bar(x,prior)
ax2.bar(x,posterior)
ax1.set_title("prior")
ax2.set_title("posterior")
ax1.set_xticks(x)
ax1.set_ylim(0,0.5)
ax2.set_ylim(0,0.5)
ax2.set_xticks(x)
plt.show()
print(tank1)
print(tank2)
for i in range(len(tank1)):plt.cla()plt.subplot(1,2,1)plt.bar(x,tank1[i])plt.xticks(x)plt.ylim(0,0.5)plt.subplot(1,2,2)plt.bar(x,tank2[i])plt.xticks(x)plt.ylim(0,0.5)plt.pause(0.5)plt.show()
#离散贝叶斯算法

import numpy as np
import random as rmstate = ["sleep","Icecream","Run"]
transitionName = [["SS","SR","SI"],["RS","RR","RI"],["IS","IR","II"]]
transitionMatrix = [[0.2,0.6,0.2],[0.1,0.6,0.3],[0.2,0.7,0.1]]
if sum(transitionMatrix[0]) + sum(transitionMatrix[1]) + sum(transitionMatrix[2]) !=3:print("Wrong")
else:print("All is gonna be okay,you should move on!!;)")def activity_forecast(days):activityToday = "Sleep"#print("Start_state:" + activityToday)activityList = [activityToday]i = 0prob = 1while i != days:if activityToday == "Sleep":change = np.random.choice(transitionName[0],replace = True,p =transitionMatrix[0])if change == "SS":prob = prob * 0.2activityList.append("Sleep")passelif change == "SR":prob = prob * 0.6activityToday = "Sleep"activityList.append("Sleep")else:prob = prob *0.2activityToday = "Icecream"activityList.append("Icecream")elif activityToday == "Run":change = np.random.choice(transitionName[1],replace = True ,p = transitionMatrix[1])if change == "RR":prob = prob * 0.5activityList.append("Run")passelif change == "RS":prob = prob *0.2activityToday = "Sleep"activityList.append("Sleep")else:prob = prob * 0.3activityToday = "Icecream"activityList.append("Icecream")elif activityToday == "Icecream":change = np.random.choice(transitionName[2],replace = True,p = transitionMatrix[2])if change == "II":prob = prob *0.1activityList.append("Icecream")passelif change == "IS":prob = prob * 0.2activityToday = "Sleep"activityList.append("Sleep")else:prob = prob * 0.7activityToday = "Run"activityList.append("Run")i += 1return activityList
list_activity = []
count = 0
for iterations in range(1,10000):list_activity.append(activity_forecast(2))
for smaller_list in list_activity:if(smaller_list[2] == "Sleep"):count += 1
percentage = (count/10000) * 100print("the probability of starting at state:'Sleep' and ending at state:'Sleep' = " + str(percentage) + "%")#马尔可夫链展示