李航老师《统计学习方法》第二版第十章课后题答案

六七～

已于 2022-12-11 17:01:13 修改

阅读量3.9k

点赞数 11

分类专栏：统计学习方法第二版文章标签：人工智能算法动态规划 python

于 2020-10-13 20:00:25 首次发布

本文链接：https://blog.csdn.net/qq_41626059/article/details/109016219

版权

统计学习方法第二版专栏收录该内容

20 篇文章 66 订阅

订阅专栏

10.1使用后向算法计算 $P(O|\lambda)$

课后题10.1
解：
话不多说，上代码，自己写的，通用版，看到网上有多个不同版本的答案，还都是手算的，很佩服这些大佬的耐心，那么多的小数，一点点的算。第一个函数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,\lambda )$

课后题10.2
解：
如果想要计算这个题的答案的话，还需要修改下上题的程序。李航老师书上有计算公式
$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 )}$
公式分析：上面公式的第一个等号成立使用了条件概率和联合概率之间的计算关系----贝叶斯公式
关键在于如何获取计算过程中的迭代向量，然后个别得到上面公式要求的，再按照公式计算就行了。

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|\lambda )=\sum_{i=1}^{N}\sum_{j=1}^{N}\alpha _{t}(i)a_{ij}b_{j}(o_{t+1})\beta _{t+1}(j)$

$t = 1, 2, . . ., T - 1$
Proof:
主要利用的是联合分布和边缘分布之间的关系
$P(O|\lambda ) = \sum_{i = 1}^{N} P(o_{1},...,o_{t-1},i_{t}=q_{i},o_{t},...,o_{T})\\ =\sum_{i=1}^{N} \sum_{j=1}^{N} P(o_{1},...,o_{t-1},i_{t}=q_{i},o_{t},i_{t+1}=q_{j},o_{t+1},...,o_{T}) \\=\sum_{i=1}^{N} \sum_{j=1}^{N} P(o_{1},...,o_{t-1},i_{t}=q_{i},o_{t})P(i_{t+1}=q_{j},o_{t+1},...,o_{T}|o_{1},...,o_{t-1},i_{t}=q_{i},o_{t}) \\=\sum_{i=1}^{N} \sum_{j=1}^{N}\alpha _{t}(i)P(i_{t+1}=q_{j},o_{t+1},...,o_{T}|i_{t}=q_{i}) \\= \sum_{i=1}^{N} \sum_{j=1}^{N}\alpha _{t}(i)P(i_{t+1}=q_{j})P(o_{t+1},...,o_{T}|i_{t}=q_{i},i_{t+1}=q_{j}) \\=\sum_{i=1}^{N} \sum_{j=1}^{N}\alpha _{t}(i)\alpha _{ij}P(o_{t+1}|i_{t+1}=q_{j})P(o_{t+2},...,o_{T}|i_{t+1}=q_{j}) \\=\sum_{i=1}^{N} \sum_{j=1}^{N}\alpha _{t}(i)\alpha _{ij}b_{j}(o_{t+1})\beta _{t+1}(j)$

10.5、比较维特比算法中变量 $\delta$ 的计算和前向算法中变量 $\alpha$ 的计算的主要区别

解：
维特比算法计算的是时刻t是给定观测序列 $o_{1},o_{2},...,o_{t}$ 计算到目前时刻t的所有状态取值的的概率最大值。
而前向算法的计算却是给定初始状态分布，计算观测序列的概率值。

六七～

关注

11
点赞
踩
20

收藏

觉得还不错? 一键收藏
6
评论
李航老师《统计学习方法》第二版第十章课后题答案

10.1解：话不多说，上代码，自己写的哦，通用版，看到网上有多个不同版本的答案，还都是手算的，也是很佩服这些大佬的耐心，那么多的小数，还一点点的算。第一个函数Bw_Recurrent(A,B,start_p,list_observation)为后向算法，名字里面的Bw表示BackWards# -*- coding: utf-8 -*-"""Created on Sun Oct 11 13:38:34 2020@author: DELL"""import numpy as npdef
复制链接

扫一扫