知乎上已經有很多的學習筆記，但讀完后總有一種這東西不是我的我理解不了的感覺，所以想試著寫一篇文章來加深一下自己的理解，也記錄下學習中的盲點。

非常推薦大家去Github看一個項目：

https://github.com/rlabbe/filterpy?github.com#下面的代碼也是完全基于上述作者的庫函數完成的，所以需要先去Github下載庫函數安裝,或者直接使用 pip install filterpy #但是作者忘記把離散的白噪聲函數放進庫里，那段函數可以在書里找到，或者改天我可以放在網盤里分享出來 #不影響下面代碼展示

這里除了有卡爾曼濾波外還有各種貝葉斯濾波的python實現，作者也寫了一本書開源在Github詳細介紹了卡爾曼，本文就算是讀書筆記吧。下面是作者發布的教材：

Jupyter Notebook Viewer?nbviewer.jupyter.org

太高估自己的效率，聽了128次大田后生仔+n首其他終于趕在沒有眼瞎前停工了。

下面代碼基本是從原作者書里粘來的，我只改了某些錯誤和可視化以及一些小的細節（動圖做的還不夠好）作者在書中試圖用追蹤狗在走廊中移動的模型來解釋離散貝葉斯濾波，由先驗概率到后驗概率再到引入反饋器這一思路詳細的介紹了貝葉斯濾波器并且實現了python的編程。同時作者一直在向我們強調學會數學編程的重要性，本人理解的不夠深刻，大家可以去翻閱原作者的書去看。

為什么我一直想實現貝葉斯濾波或者卡爾曼濾波，因為它是一種動態濾波器，它可以用來建立動態的反演模型，這和我以前接觸的反演很不一樣，我曾經做過利用光學遙感建立泥沙反演模型和利用兩層BP神經網絡對地震波建立反演模型找地下礦層的兩個小東西，以上兩個反演模型都是靜態的過程，而卡爾曼濾波最反智的則是在濾波中不斷更新模型參數（加入了后饋過程也就是K卡爾曼增益矩陣），這讓我困惑了特別久。

如果大家也對這種濾波器感興趣建議先學些統計知識和最小二乘法，卡爾曼濾波和最小二乘法之間有著很緊密地聯系，并且最小二乘法方程可以很好的幫助你理解最小均方誤差，這是實現卡爾曼濾波的重要假設。

下面也會放一段簡單馬爾可夫鏈的python小程序，如果你和我一樣用Sublime ，多試幾次Ctrl + B，這是基于隨機數做的，但是你會發現在多次迭代后概率會收斂到一個定值，還記得統計課上講的大數定律嘛，在樣本空間足夠大時，依概率收斂到某個數，大約就是這個意思。所以在隨機天氣生成器中需要用蒙托卡羅思想來添加擾動來破壞這種收斂性？，我上一篇文章中也寫道過，但是我想錯了，并不是對轉移概率矩陣添加擾動，而是對初始狀態添加擾動（最近聽了幾個MCMC方法的讀書匯報，發現了這個秘密）

有關于今天統計課上講到的氣候序列均一化我也想說點兒給某位不來上課的同學，我認為卡爾曼濾波器可以很好的應用到均一化過程中，具體方法我不說，啦啦啦。

#離散貝葉斯濾波器的實現展示 #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) + "%")#馬爾可夫鏈展示

能看到這兒說明你對我的文章很有興趣，那么有任何問題可以通過任何方式問我。

日常吐槽：

最近效率真是低的可怕，可我睡不著啊。。。

休息不好吃又吃不下要掉體重了，為了回家不讓老媽發現（不然要被瘋狂投食）牛奶得換全脂來保證熱量攝入了。

來吧，都看到這兒了陪我聽首歌啊：

キミがいれば(世紀末ヴァージョン)?music.163.com

柯姓男孩永不認輸

總結

以上是生活随笔為你收集整理的python 排序统计滤波器_马尔可夫链+贝叶斯滤波器的Python展示的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。