10.1使用后向演算法計算 P ( O ∣ λ ) P(O|\lambda) P(O∣λ)

解:
話不多說,上代碼,自己寫的,通用版,看到網上有多個不同版本的答案,還都是手算的,很佩服這些大佬的耐心,那么多的小數,一點點的算,第一個函式Bw_Recurrent(A,B,start_p,list_observation)為后向演算法,名字里面的Bw表示BackWards
# -*- coding: utf-8 -*-
import numpy as np
def Bw_Recurrent(A,B,start_p,list_observation):
'''
Parameters
----------
A : np.float32 matrix
是隱馬模型的狀態轉移矩陣.
B : np.float32 matrix
觀測概率矩陣.
start_p : np.float32 matrix
初始狀態概率分布.
list_observation:list
用于計算特定序列的概率
Returns
-------
float32 .
是P(O|lamda)
'''
#后向的遞回計算需要初始化一個全一的大小為N*1的矩陣
N=np.shape(A)[0]
p=np.ones((N,1),dtype=np.float32)
T=len(list_observation)#觀測序列的長度
for i in range(T-1):
#將觀測矩陣里面的第list_observation[T-i-1]列取出來
v=np.transpose(np.array([B[:,list_observation[T-i-1]]],dtype=np.float32))
p=np.matmul(A,v*p)#這行代碼既有矩陣乘法,也有矩陣點乘
o1=np.transpose(np.array([B[:,list_observation[0]]],dtype=np.float32))
return np.sum(start_p*o1*p)#注意里面的矩陣乘法是點乘操作,也就是將對應位置的元素相乘
def Fw_Recurrent(A,B,start_p,list_observation):
"""
Parameters
----------
A : np.float32 matrix
是隱馬模型的狀態轉移矩陣.
B : np.float32 matrix
觀測概率矩陣.
start_p : np.float32 matrix
初始狀態概率分布.
list_observation:list
用于計算特定序列的概率
Returns
-------
float32 .
是P(O|lamda)
"""
o1=np.transpose(np.array([B[:,list_observation[0]]],dtype=np.float32))
p=start_p*o1
T=len(list_observation)#觀測序列的長度
for i in range(T-1):
p=np.matmul(np.transpose(A),p)*np.transpose(np.array([B[:,list_observation[i+1]]],dtype=np.float32))
#要注意上面這行代碼里面的矩陣乘法和矩陣點乘
print(p,'\n')
return np.sum(p)
if __name__=='__main__':
A=np.array([[0.5,0.2,0.3],
[0.3,0.5,0.2],
[0.2,0.3,0.5]],dtype=np.float32)
B=np.array([[0.5,0.5],
[0.4,0.6],
[0.7,0.3]],dtype=np.float32)
start_p=np.array([[0.2],[0.4],[0.4]],dtype=np.float32)
list_observation=[0,1,0,1]
s_fw=Fw_Recurrent(A,B,start_p,list_observation)
s_bw=Bw_Recurrent(A,B,start_p,list_observation)
print('前向演算法的結果是: ',s_fw)
print('后向演算法的結果是: ',s_bw)
運行結果是:
前向演算法的結果是: 0.060090814
后向演算法的結果是: 0.060090803
由上圖可以看到,這里輸出的結果有十分微小的差距,這是因為計算機的截斷誤差導致的
10.2 使用前后向演算法計算 P ( i 4 = q 3 ∣ O , λ ) P(i_{4}=q_{3}|O,\lambda ) P(i4?=q3?∣O,λ)

解:
如果想要計算這個題的答案的話,還需要修改下上題的程式,李航老師書上有計算公式
P
(
i
t
=
q
i
∣
O
,
λ
)
=
P
(
i
t
=
q
i
,
O
∣
λ
)
P
(
O
∣
λ
)
=
α
t
(
i
)
?
β
t
(
i
)
P
(
O
∣
λ
)
P(i_{t}=q_{i}|O,\lambda )=\frac{P(i_{t}=q_{i},O|\lambda )}{P(O|\lambda )}=\frac{\alpha _{t}(i)*\beta _{t}(i)}{P(O|\lambda )}
P(it?=qi?∣O,λ)=P(O∣λ)P(it?=qi?,O∣λ)?=P(O∣λ)αt?(i)?βt?(i)?
公式分析:上面公式的第一個等號成立使用了條件概率和聯合概率之間的計算關系----貝葉斯公式
關鍵在于如何獲取計算程序中的迭代向量,然后個別得到上面公式要求的,再按照公式計算就行了,
import numpy as np
def Bw_Recurrent(A,B,start_p,list_observation):
'''
Parameters
----------
A : np.float32 matrix
是隱馬模型的狀態轉移矩陣.
B : np.float32 matrix
觀測概率矩陣.
start_p : np.float32 matrix
初始狀態概率分布.
list_observation:list
用于計算特定序列的概率
Returns
-------
float32 .
是P(O|lamda)
'''
P=[]
#后向的遞回計算需要初始化一個全一的大小為N*1的矩陣
N=np.shape(A)[0]
p=np.ones((N,1),dtype=np.float32)
P.append(p)
T=len(list_observation)#觀測序列的長度
for i in range(T-1):
#將觀測矩陣里面的第list_observation[T-i-1]列取出來
v=np.transpose(np.array([B[:,list_observation[T-i-1]]],dtype=np.float32))
p=np.matmul(A,v*p)#這行代碼既有矩陣乘法,也有矩陣點乘
P.append(p)
o1=np.transpose(np.array([B[:,list_observation[0]]],dtype=np.float32))
return np.sum(start_p*o1*p),P#注意里面的矩陣乘法是點乘操作,也就是將對應位置的元素相乘
def Fw_Recurrent(A,B,start_p,list_observation):
"""
Parameters
----------
A : np.float32 matrix
是隱馬模型的狀態轉移矩陣.
B : np.float32 matrix
觀測概率矩陣.
start_p : np.float32 matrix
初始狀態概率分布.
list_observation:list
用于計算特定序列的概率
Returns
-------
float32 .
是P(O|lamda)
"""
P=[]
o1=np.transpose(np.array([B[:,list_observation[0]]],dtype=np.float32))
p=start_p*o1
P.append(p)
T=len(list_observation)#觀測序列的長度
for i in range(T-1):
p=np.matmul(np.transpose(A),p)*np.transpose(np.array([B[:,list_observation[i+1]]],dtype=np.float32))
#要注意上面這行代碼里面的矩陣乘法和矩陣點乘
P.append(p)
return np.sum(p),P
if __name__=='__main__':
A=np.array([[0.5,0.1,0.4],
[0.3,0.5,0.2],
[0.2,0.2,0.6]],dtype=np.float32)
B=np.array([[0.5,0.5],
[0.4,0.6],
[0.7,0.3]],dtype=np.float32)
start_p=np.array([[0.2],[0.3],[0.5]],dtype=np.float32)
list_observation=[0,1,0,0,1,0,1,1]
s_fw,P_fw=Fw_Recurrent(A,B,start_p,list_observation)
s_bw,P_bw=Bw_Recurrent(A,B,start_p,list_observation)
print('前向演算法的結果是: ',s_fw)
print('后向演算法的結果是: ',s_bw)
print('result is :',P_fw[3][2]*P_bw[4][2]/s_bw)
計算結果是:
前向演算法的結果是: 0.00347671
后向演算法的結果是: 0.0034767103
result is : [0.5369518]
10.3 在習題10.1中,試用維特比(Viterbi)演算法求最優的路徑
解:
看到維特比演算法在很多地方都有應用,但是在不同的地方應用的方式是不同的,如果給了柵欄網圖,并且再給每一條邊都附上權值,我相信很多都會計算,但是,不同的應用場景中,如何計算這個權值是個關鍵,根據書上列好的演算法無論是寫程式還是計算起來都很簡單,但是,我覺得關鍵是如何做到活學活用,在不同的場景里面應用自如,好藍呀!!!
一個演算法學了差不多一個下午,但是感徑訓是差點感性認識,就這樣吧,上代碼
# -*- coding: utf-8 -*-
"""
Created on Mon Oct 12 15:52:10 2020
@author: DELL
"""
import numpy as np
def viterbi_decode(A,B,start_p,list_observation):
'''
Parameters
----------
A : 2D matrix float
狀態轉移矩陣.
B : 2D matrix float
觀測概率矩陣.
start_p : float matrix
初始狀態概率分布.
observation_list : list
觀測序列.
Returns
-------
最優路徑.
'''
best_path=[]
T=len(list_observation)
#初始化
N=A.shape[0]#狀態個數
O1=np.transpose(np.array([B[:,list_observation[0]]],dtype=np.float32))
delta=start_p*O1
print(delta,'\n')
all_psi=[]
psi=np.zeros([N,1],dtype=np.float32)
for i in range(T-1):
psi=np.argmax(delta*A,axis=0)
delta=np.transpose(np.array([np.max(delta*A,axis=0)]))*np.transpose(np.array([B[:,list_observation[i+1]]],dtype=np.float32))
all_psi.append(psi)
final_score=np.max(delta)
best_path.append(np.argmax(delta))
for i in range(T-1):
best_path.insert(0,all_psi[T-2-i][best_path[0]])
return final_score,best_path
if __name__=='__main__':
A=np.array([[0.5,0.2,0.3],
[0.3,0.5,0.2],
[0.2,0.3,0.5]],dtype=np.float32)
B=np.array([[0.5,0.5],
[0.4,0.6],
[0.7,0.3]],dtype=np.float32)
start_p=np.array([[0.2],[0.4],[0.4]],dtype=np.float32)
list_observation=[0,1,0,1]
final_score,best_path=viterbi_decode(A,B,start_p,list_observation)
print('最優路徑的概率是:',final_score)
print('最優路徑為:',best_path)
運行結果是:
最優路徑的概率是: 0.0030240004
最優路徑為: [2, 1, 1, 1]
分析:因為我的狀態代碼是從0開始的,因此對應起來的正確路徑就是[3,2,2,2]
書上的例子運行結果是:
最優路徑的概率是: 0.014700001
最優路徑為: [2, 2, 2]
對應起來就是[3,3,3]
10.4、試用前向概率和后向概率推導 P ( O ∣ λ ) = ∑ i = 1 N ∑ j = 1 N α t ( i ) a i j b j ( o t + 1 ) β t + 1 ( j ) P(O|\lambda )=\sum_{i=1}^{N}\sum_{j=1}^{N}\alpha _{t}(i)a_{ij}b_{j}(o_{t+1})\beta _{t+1}(j) P(O∣λ)=i=1∑N?j=1∑N?αt?(i)aij?bj?(ot+1?)βt+1?(j)
t
=
1
,
2
,
.
.
.
,
T
?
1
t=1,2,...,T-1
t=1,2,...,T?1
Proof:
根據后向演算法的遞推公式計算可得:
β
t
(
i
)
=
∑
j
=
1
N
a
i
j
b
j
(
o
t
+
1
)
β
t
+
1
(
j
)
i
=
1
,
2
,
.
.
.
,
N
(1)
\beta _{t}(i)=\sum_{j=1}^{N}a_{ij}b_{j}(o_{t+1})\beta _{t+1}(j)\quad i=1,2,...,N\tag{1}
βt?(i)=j=1∑N?aij?bj?(ot+1?)βt+1?(j)i=1,2,...,N(1)
又根據李航老師書上的公式可得:
P
(
O
∣
λ
)
=
∑
j
=
1
N
α
t
(
j
)
β
t
(
j
)
(2)
P(O|\lambda )=\sum_{j=1}^{N}\alpha _{t}(j)\beta _{t}(j) \tag{2}
P(O∣λ)=j=1∑N?αt?(j)βt?(j)(2)
就算不使用李航老師書上的公式,也可以很簡單的推出公式(2)
下面推導一下:
P
(
O
∣
λ
)
=
∑
i
=
1
N
∑
j
=
1
N
α
t
(
i
)
a
i
j
b
j
(
o
t
+
1
)
β
t
+
1
(
j
)
=
∑
i
=
1
N
α
t
(
i
)
∑
j
=
1
N
a
i
j
b
j
(
o
t
+
1
)
β
t
+
1
(
j
)
P(O|\lambda )=\sum_{i=1}^{N}\sum_{j=1}^{N}\alpha _{t}(i)a_{ij}b_{j}(o_{t+1})\beta _{t+1}(j)=\sum_{i=1}^{N}\alpha _{t}(i)\sum_{j=1}^{N}a_{ij}b_{j}(o_{t+1})\beta _{t+1}(j)
P(O∣λ)=i=1∑N?j=1∑N?αt?(i)aij?bj?(ot+1?)βt+1?(j)=i=1∑N?αt?(i)j=1∑N?aij?bj?(ot+1?)βt+1?(j)
=
∑
i
=
1
N
α
t
(
i
)
β
t
(
i
)
=\sum_{i=1}^{N}\alpha _{t}(i)\beta _{t}(i)
=i=1∑N?αt?(i)βt?(i)
下面給出詳細的證明,(2)式的推導需要用到聯合分布與邊際分布之間的關系,還有最重要的就是一個假設:觀測獨立性假設,也就是任意時刻的觀測值只依賴于該時刻馬爾可夫鏈的狀態值,與其他狀態和觀測值無關,下面的公式(3)成立,
P
(
o
t
∣
i
T
,
o
T
,
i
T
?
1
,
o
T
?
1
,
.
.
.
,
i
1
,
o
1
)
=
p
(
o
t
∣
i
t
)
(3)
P(o_{t}|i_{T},o_{T},i_{T-1},o_{T-1},...,i_{1},o_{1})=p(o_{t}|i_{t})\tag{3}
P(ot?∣iT?,oT?,iT?1?,oT?1?,...,i1?,o1?)=p(ot?∣it?)(3)
根據聯合分布于邊際分布之間的關系,那么就有:
P
(
O
∣
λ
)
=
∑
i
=
1
N
P
(
O
,
i
t
=
q
i
∣
λ
)
(4)
P(O|\lambda )=\sum_{i=1}^{N}P(O,i_{t}=q_{i}|\lambda )\tag{4}
P(O∣λ)=i=1∑N?P(O,it?=qi?∣λ)(4)
上式的意義也就是對所有的在
i
t
i_{t}
it?時刻的狀態求和得到邊際分布
P
(
O
∣
λ
)
P(O|\lambda )
P(O∣λ)
再根據前向演算法后向演算法的定義,所以有:
P
(
O
∣
λ
)
=
∑
i
=
1
N
P
(
O
,
i
t
=
q
i
∣
λ
)
=
∑
i
=
1
N
P
(
o
1
,
o
2
,
.
.
.
,
o
t
,
.
.
.
,
o
T
,
i
t
=
q
t
∣
λ
)
P(O|\lambda)=\sum_{i=1}^{N}P(O,i_{t}=q_{i}|\lambda )=\sum_{i=1}^{N}P(o_{1},o_{2},...,o_{t},...,o_{T},i_{t}=q_{t}|\lambda )
P(O∣λ)=i=1∑N?P(O,it?=qi?∣λ)=i=1∑N?P(o1?,o2?,...,ot?,...,oT?,it?=qt?∣λ)
=
∑
i
=
1
N
P
(
o
1
,
o
2
,
.
.
.
,
o
t
,
i
t
=
q
t
∣
λ
)
?
P
(
o
t
+
1
,
.
.
.
,
o
T
∣
o
1
,
o
2
,
.
.
.
,
o
t
,
i
t
=
q
t
,
λ
)
=\sum_{i=1}^{N}P(o_{1},o_{2},...,o_{t},i_{t}=q_{t}|\lambda )*P(o_{t+1},...,o_{T}|o_{1},o_{2},...,o_{t},i_{t}=q_{t},\lambda)
=i=1∑N?P(o1?,o2?,...,ot?,it?=qt?∣λ)?P(ot+1?,...,oT?∣o1?,o2?,...,ot?,it?=qt?,λ)
=
∑
i
=
1
N
P
(
o
1
,
o
2
,
.
.
.
,
o
t
,
i
t
=
q
t
∣
λ
)
?
P
(
o
t
+
1
,
.
.
.
,
o
T
∣
i
t
=
q
t
,
λ
)
=\sum_{i=1}^{N}P(o_{1},o_{2},...,o_{t},i_{t}=q_{t}|\lambda )*P(o_{t+1},...,o_{T}|i_{t}=q_{t},\lambda)
=i=1∑N?P(o1?,o2?,...,ot?,it?=qt?∣λ)?P(ot+1?,...,oT?∣it?=qt?,λ)
=
∑
j
=
1
N
α
t
(
j
)
β
t
(
j
)
t
=
1
,
2
,
.
.
.
,
T
?
1
=\sum_{j=1}^{N}\alpha _{t}(j)\beta _{t}(j) \quad t=1,2,...,T-1
=j=1∑N?αt?(j)βt?(j)t=1,2,...,T?1
得證,
獨立性的假設是一個比較強的條件,其實作實往往很復雜,這里只是簡化了條件,不然估計很多理論難以實作到實際中去,
10.5、比較維特比演算法中變數 δ \delta δ的計算和前向演算法中變數 α \alpha α的計算的主要區別
解:
維特比演算法計算的是時刻t是給定觀測序列
o
1
,
o
2
,
.
.
.
,
o
t
o_{1},o_{2},...,o_{t}
o1?,o2?,...,ot?計算到目前時刻t的所有狀態取值的的概率最大值,
而前向演算法的計算卻是給定初始狀態分布,計算觀測序列的概率值,
轉載請註明出處,本文鏈接:https://www.uj5u.com/qita/172389.html
標籤:其他
下一篇:python 使用攝像頭監測心率
