《数据挖掘》第六课笔记


一、蒙特卡罗思想之接受-拒绝采样

  在计算积分值时,蒙特拉洛的主要内容就是从假设分布中抽取n次,估计积分值。对于概率分布不是常见的分布,一个可行的办法是采用接受-拒绝采样来得到该分布的样本。既然 p ( x ) p(x) p(x)太复杂在程序中没法直接采样,那么我设定一个程序可采样的分布 q ( x ) q(x) q(x) 比如正态分布,然后按照一定的方法拒绝某些样本,以达到接近 p ( x ) p(x) p(x) 分布的目的。
在这里插入图片描述
  具体采用过程如下,设定一个方便采样的常用概率分布函数 q ( x ) q(x) q(x),以及一个常量 k k k,使得 p ( x ) p(x) p(x)总在 k q ( x ) kq(x) kq(x) 的下方,如上图
  首先,采样得到 q ( x ) q(x) q(x)的一个样本 z 0 {z_0} z0,然后,从均匀分布 ( 0 , k q ( z 0 ) ) (0,kq({z_0})) (0,kq(z0))中采样得到一个值 u u u。如果 u u u落在了上图中的灰色区域,则拒绝这次抽样,否则接受这个样本 z 0 {z_0} z0。重复以上过程得到 n n n个接受的样本 z 0 , z 1 , … z n − 1 {z_0},{z_1}, \ldots {z_{n - 1}} z0,z1,zn1,则最后的蒙特卡罗求解结果为: p ( z i ) k q ( z i ) {\frac{{p({z_i})}}{{kq({z_i})}}} kq(zi)p(zi)  整个过程中,我们通过一系列的接受拒绝决策来达到用 q ( x ) q(x) q(x)模拟 p ( x ) p(x) p(x)概率分布的目的。


二、马尔可夫链

1.马尔可夫链定义

  随机过程 { X n , n = 0 , 1 , 2 , … L } \{{X_{\rm{n}}},n = 0,1,2, \ldots L\} {Xn,n=0,1,2,L}称为Markov链,若它只取有限值或可列个值(称为过程的状态,记为0,1,2,…),{0,1,2…}或者其子集记为S,称为过程的状态空间,对任意 n n n ≥ \ge 0 0 0,及状态 i , j , i 0 , i 1 , L , i n − 1 i,j,{i_0},{i_1},L,{i_{n - 1}} i,j,i0,i1,L,in1,有(马氏性) P ( X n + 1 = j ∣ X 0 = i 0 , X 1 = i 1 , L , X n − 1 = i n − 1 , X n = i ) = P ( X n + 1 = j ∣ X n = i ) \color{red}P({X_{n + 1}} = j|{X_0} = {i_0},{X_1} = {i_1},L,{X_{n - 1}} = {i_{n - 1}},{X_n} = i) = P({X_{n + 1}} = j|{X_n} = i) P(Xn+1=jX0=i0,X1=i1,L,Xn1=in1,Xn=i)=P(Xn+1=jXn=i)  即马尔可夫链 { X n , n ≥ 0 } \{ {X_{\rm{n}}},n \ge 0\} {Xn,n0}的有限维分布完全由初始分布 P { X 0 = i 0 } P\{ {X_0} = {i_0}\} P{X0=i0}条件概率(一步转移概率) P ( X n = j ∣ X n − 1 = i ) P({X_{n }} = j|{X_{n-1}} = i) P(Xn=jXn1=i)确定

2.马尔可夫链状态的分类和性质

2.1 状态的分类

  定义1:若存在 n n n ≥ \ge 0 0 0使得 P i j ( n ) P_{ij}^{(n)} Pij(n) > \gt > 0 0 0,则称从状态 i i i可达状态 j j j,记 i → j i \to j ij。如果 i → j i \to j ij j → i j \to i ji,则称 i i i j j j互通,记 i ↔ j \color{red}i \leftrightarrow j ij
  将任何两个互通的状态归为一类:(1)同一类的状态之间都是互通的.(2)任何一个状态不能同时属于两个不同的类。
  定义2:若Markov链只存在一个类,就称它是不可约的,否则称为可约的

2.2 状态的性质

  定义3:若集合 { n : n ≥ 1 , P i i ( n ) > 0 } ≠ ∅ \{ n:n \ge 1,P_{ii}^{(n)} > 0\} \ne \emptyset {n:n1,Pii(n)>0}=,则称它的最大公约数 d = d ( j ) d=d(j) d=d(j)为状态 i i i的周期,若 { d = 1 , 称 i 为非周期的 d > 1 , 称 i 为周期的 \color{red}\{ _{d = 1,称i为非周期的}^{d > 1,称i为周期的} {d=1,i为非周期的d>1,i为周期的特.若上述集合为空集,则称$i$的周期为无穷大
  定理:若状态 i , j i,j i,j同属一类,则 d ( i ) = d ( j ) d(i)=d(j) d(i)=d(j)

2.3 平稳分布(不变分布)与极限分布

  平稳分布:设马氏链 { X n , n ≥ 0 } \{ {X_{\rm{n}}},n \ge 0\} {Xn,n0}有转移概率矩阵 P = ( P i j ) P{\rm{ = (}}{{\rm{P}}_{{\rm{ij}}}}{\rm{)}} P=(Pij),若存在一个概率分布 { π i , i ∈ S } \{ {\pi _i},i \in S\} {πi,iS},满足 π j = ∑ i ∈ S π i P i j {\pi _j} = \sum\limits_{i \in S} {{\pi _i}{P_{ij}}} πj=iSπiPij,则称 { π i , i ∈ S } \color{red}\{ {\pi _i},i \in S\} {πi,iS}为该链的平稳分布(不变分布)
  分析:
   令 π = ( π 1 , π 2 , … ) \pi = ({\pi _1},{\pi _2}, \ldots ) π=(π1,π2,),则上式为 π = π P \pi = \pi P π=πP
   (1).由 π = π P \pi = \pi P π=πP可知, π ( I − P ) = 0 \pi (I - P) = 0 π(IP)=0,故 I I I是矩阵 P P P的左特征值,平稳分布 π \pi π P P P的左特征向量
   (2).两边同乘 P P P,得 π = π P = π P 2 = … = π P n \pi = \pi P = \pi {P^2} = \ldots = \pi {P^n} π=πP=πP2==πPn

  极限分布:称Markov链是遍历的,如果所有状态相通且均是周期为1的正常返状态,对于遍历的Markov链,极限 lim ⁡ n → ∞ P i j ( n ) = π j = 1 u j \color{red}\mathop {\lim }\limits_{n \to \infty } P_{ij}^{(n)} = {\pi _j} = \frac{1}{{{u_j}}} nlimPij(n)=πj=uj1 j ∈ S j \in S jS称为Markov链的极限分布

3.基于马尔可夫链采样

  如果我们得到了某个平稳分布所对应的马尔科夫链状态转移矩阵 P P P,我们就很容易采用出这个平稳分布的样本集。
  假设我们任意初始的概率分布是 π 0 ( x ) \pi_0(x) π0(x),经过第一轮马尔科夫链状态转移后的概率分布是 π 1 ( x ) \pi_1(x) π1(x),第 i i i轮的概率分布 π i ( x ) \pi_i(x) πi(x)。 假设经过 n n n轮后马尔科夫链收敛到我们的平稳分布 π ( x ) \pi(x) π(x),即: π n ( x ) = π n + 1 ( x ) = π n + 2 ( x ) = … = π ( x ) {\pi _n}(x) = {\pi _{n + 1}}(x) = {\pi _{n + 2}}(x) = \ldots = \pi (x) πn(x)=πn+1(x)=πn+2(x)==π(x)  对于每个分布 π i ( x ) \pi_i(x) πi(x),我们有: π i ( x ) = π i − 1 ( x ) P = π i − 2 ( x ) P 2 = … = π 0 ( x ) P i {\pi _i}(x) = {\pi _{i - 1}}(x)P = {\pi _{i - 2}}(x){P^2} = \ldots = {\pi _0}(x){P^i} πi(x)=πi1(x)P=πi2(x)P2==π0(x)Pi  现在我们可以开始采样了,首先,基于初始任意简单概率分布比如高斯分布 π 0 ( x ) \pi_0(x) π0(x)采样得到状态值 x 0 x_0 x0,基于条件概率分布 P ( x ∣ x 0 ) P(x|x_0) P(xx0)采样状态值 x 1 x_1 x1 ,一直进行下去,当状态转移进行到一定的次数时,比如到 n n n次时,我们认为此时的采样集 ( x n , x n + 1 , x n + 2 , . . . ) (x_n,x_{n+1},x_{n+2},...) (xn,xn+1,xn+2,...)即是符合我们的平稳分布的对应样本集,可以用来做蒙特卡罗模拟求和了。

  总结下基于马尔科夫链的采样过程:
  (1)输入马尔科夫链状态转移矩阵 P P P,设定状态转移次数阈值 n 1 n_1 n1 ,需要的样本个数 n 2 n_2 n2
  (2)从任意简单概率分布采样得到初始状态值 x 0 x_0 x0
  (3) f o r   t = 0   t o   n 1 + n 2 − 1 for\ t=0\ to\ n_1+n_2−1 for t=0 to n1+n21: 从条件概率分布 P ( x ∣ x t ) P(x|x_t) P(xxt)中采样得到样本 x t + 1 x_{t+1} xt+1
  (4)样本集 ( x n 1 , x n 1 + 1 , . . . , x n 1 + n 2 − 1 ) \color{red}(x_{n_1}, x_{n_1+1},..., x_{n_1+n_2-1}) (xn1,xn1+1,...,xn1+n21)即为我们需要的平稳分布对应的样本集。

4.一阶马尔科夫链模型(牛市熊市股票例子)

链接: 《数据挖掘》第二次实验


三、隐马尔可夫链基本定义

1.隐马氏模型的组成和基本假设

  组成:
  • 初始概率分布
  • 状态转移概率分布
  • 观测概率分布
  • Q:所有可能状态的集合
  • V:所有可能观测的集合
  • I: 长度为T的状态序列
  • O:对应的观测序列
  • A:状态转移概率矩阵
  • B:观测概率矩阵
  • π \pi π:初始状态概率向量
  基本假设:
  (1)齐次马尔科夫性假设,隐马尔可分链t的状态只和t-1状态有关。
  (2)观测独立性假设,观测只和当前时刻状态有关。

2.隐马氏模型的3个基本问题

(1)概率计算问题
  给定: λ = ( A , B , π ) \lambda = (A,B,\pi ) λ=(A,B,π)   O = ( o 1 , o 2 , ⋯   , o T ) O = ({o_1},{o_2}, \cdots ,{o_T}) O=(o1,o2,,oT)
  计算: P ( O ∣ λ ) P(O|\lambda ) P(Oλ)
(2)学习问题
  已知: O = ( o 1 , o 2 , ⋯   , o T ) O = ({o_1},{o_2}, \cdots ,{o_T}) O=(o1,o2,,oT)
  估计: λ = ( A , B , π ) \lambda = (A,B,\pi ) λ=(A,B,π),使 P ( O ∣ λ ) P(O|\lambda ) P(Oλ)最大
(3)预测问题(解码)
  已知: λ = ( A , B , π ) \lambda = (A,B,\pi ) λ=(A,B,π)   O = ( o 1 , o 2 , ⋯   , o T ) O = ({o_1},{o_2}, \cdots ,{o_T}) O=(o1,o2,,oT)
  求:使 P ( O ∣ λ ) P(O|\lambda ) P(Oλ)最大的状态序列 I = ( i 1 , i 2 , ⋯   , i T ) I = ({i_1},{i_2}, \cdots ,{i_T}) I=(i1,i2,,iT)

3.隐马氏模型的应用

  • 人脸识别
  • 语音识别
  • 入侵检测
  • 情报领域
  • 手写体识别
  • 通信领域的解码器


四、隐马尔可夫链的概率计算问题求解方法

1.暴力求解

  1.一个简单问题
  知道骰子有几种,每种骰子是什么,每次掷的都是什么骰子,根据掷骰子掷出的结果,求产生这个结果的概率。
在这里插入图片描述
P = P ( D 6 ) ∗ P ( D 6 → 1 ) ∗ P ( D 6 → D 8 ) ∗ P ( D 8 → 6 ) ∗ P ( D 8 → D 8 ) ∗ P ( D 8 → 3 ) P = P(D6)*P(D6 \to 1)*P(D6 \to D8)*P(D8 \to 6)*P(D8 \to D8)*P(D8 \to 3) P=P(D6)P(D61)P(D6D8)P(D86)P(D8D8)P(D83)
= 1 3 ∗ 1 6 ∗ 1 3 ∗ 1 8 ∗ 1 3 ∗ 1 8 = \frac{1}{3}*\frac{1}{6}*\frac{1}{3}*\frac{1}{8}*\frac{1}{3}*\frac{1}{8} =316131813181

  2.谁动了我的骰子?
  比如说你怀疑自己的六面骰被赌场动过手脚了,有可能被换成另一种六面骰,这种六面骰掷出来是1的概率更大,是1/2,掷出来是2,3,4,5,6的概率是1/10。你怎么办么?
  答案很简单,算一算正常的三个骰子掷出一段序列的概率,再算一算不正常的六面骰和另外两个正常骰子掷出这段序列的概率。如果前者比后者小,你就要小心了。比如说掷骰子的结果是:
在这里插入图片描述
  简单而暴力的方法就是把穷举所有的骰子序列,还是计算每个骰子序列对应的概率相加,得到的总概率就是我们要求的结果。
  首先,如果我们只掷一次骰子:
在这里插入图片描述
  看到结果为1。产生这个结果的总概率可以按照如下计算,总概率为0.18:
在这里插入图片描述
  把这个情况拓展,我们掷两次骰子在这里插入图片描述
  看到结果为1,6。产生这个结果的总概率可以按照如下计算,总概率为0.05:
在这里插入图片描述
  继续拓展,我们掷三次骰子:在这里插入图片描述
  看到结果为1,6,3。产生这个结果的总概率可以按照如下计算,总概率为0.03:在这里插入图片描述

  暴力求解算法的时间复杂度是 O ( T N T ) O(T{N^T}) O(TNT)

2.前向算法

  仅做一个基本介绍,前向算法的时间复杂度是 O ( N 2 T ) O({N^2}T) O(N2T)
  前向概率定义:给定隐马尔可夫模型
定义到时刻 t t t部分观测序列: o 1 , o 2 , ⋯   , o i {o_1},{o_2}, \cdots ,{o_i} o1,o2,,oi,且状态为 q i {q_i} qi的概率为前向概率,记作 α t ( i ) = P ( o 1 , o 2 , ⋯   , o t , i t = q i ∣ λ ) {\alpha _t}(i) = P({o_1},{o_2}, \cdots ,{o_t},{i_t} = {q_i}|\lambda ) αt(i)=P(o1,o2,,ot,it=qiλ)
  算法(观测序列概率的前向算法)
  输入:隐马尔可夫模型 λ \lambda λ,观测序列 O O O
  输出:观测序列概率 P ( O ∣ λ ) P(O|\lambda ) P(Oλ)
  初值: α t ( i ) = π i b i ( o 1 ) {\alpha _t}(i) = {\pi _i}{b_{\rm{i}}}({o_1}) αt(i)=πibi(o1) i = 1 , 2 , . . . , N i=1,2,...,N i=1,2,...,N
  递推: α t + 1 ( i ) = [ ∑ i = 1 N α t ( j ) α j i ] b i ( o t + 1 ) {\alpha _{t + 1}}(i) = \left[ {\sum\limits_{i = 1}^N {{\alpha _t}(j){\alpha _{ji}}} } \right]{b_i}({o_{t + 1}}) αt+1(i)=[i=1Nαt(j)αji]bi(ot+1) i = 1 , 2 , . . . , N i=1,2,...,N i=1,2,...,N
  终止: P ( O ∣ λ ) = ∑ i = 1 N α T ( i ) P(O|\lambda ) = \sum\limits_{i = 1}^N {{\alpha _T}(i)} P(Oλ)=i=1NαT(i)


五、隐马尔可夫链的预测问题求解办法之viterbi维特比算法

1.算法(维特比算法)

  输入:模型 λ = ( A , B , π ) \lambda = (A,B,\pi ) λ=(A,B,π)和观测 O = ( o 1 , o 2 , ⋯   , o T ) O = ({o_1},{o_2}, \cdots ,{o_T}) O=(o1,o2,,oT)
  输出:最优路径 I ∗ = ( i 1 ∗ , i 2 ∗ , ⋯   , i T ∗ ) {I^*} = (i_1^*,i_2^*, \cdots ,i_T^*) I=(i1,i2,,iT)
  (1)初始化 δ 1 ( i ) = π i b i ( o 1 ) , i = 1 , 2 , . . . N {\delta _1}(i) = {\pi _i}{b_i}({o_1}),i=1,2,...N δ1(i)=πibi(o1)i=1,2,...N ϕ 1 ( i ) = 0 , i = 1 , 2 , . . . N {\phi _1}(i) = 0,i=1,2,...N ϕ1(i)=0i=1,2,...N  (2)递推.对 t = 2 , 3 , . . . , T t=2,3,...,T t=2,3,...,T δ t ( i ) = max ⁡ 1 ≤ j ≤ N [ δ t − 1 ( j ) a j i ] b i ( o t ) , i = 1 , 2 , . . . N {\delta _t}(i) = \mathop {\max }\limits_{1 \le j \le N} [{\delta _{t - 1}}(j){a_{ji}}]{b_i}({o_t}),i=1,2,...N δt(i)=1jNmax[δt1(j)aji]bi(ot)i=1,2,...N ϕ t ( i ) = arg ⁡ max ⁡ 1 ≤ j ≤ N [ δ t − 1 ( j ) a j i ] , i = 1 , 2 , . . . N {\phi _t}(i) = \arg \mathop {\max }\limits_{1 \le j \le N} [{\delta _{t - 1}}(j){a_{ji}}],i=1,2,...N ϕt(i)=arg1jNmax[δt1(j)aji]i=1,2,...N

2.盒球模型

在这里插入图片描述

3.东京天气预测问题

  一个东京的朋友每天根据天气{下雨,天晴}决定当天的活动{公园散步,购物,清理房间}中的一种,我每天只能在twitter上看到她发的动态“啊,我前天公园散步、昨天购物、今天清理房间了!”,那么我可以根据她发的推特推断东京这三天的天气。
  在这个例子里,显状态是活动,隐状态是天气。

3.1手写计算

在这里插入图片描述

3.2 code

import numpy as np

class HMM():
    def __init__(self):
        # 定义天气状态和对应的索引
        self.states = {'下雨': 0, '天晴': 1}
        # 定义观测状态和对应的索引
        self.observations = {'公园散步': 0, '购物': 1, '清理房间': 2}
        # 定义初始状态概率
        self.pai = np.array([0.6, 0.4])
        # 定义状态转移概率
        self.A = np.array([[0.7, 0.3], [0.4, 0.6]])
        # 定义观测概率
        self.B = np.array([[0.1, 0.4, 0.5], [0.6, 0.3, 0.1]])
        # 观测序列
        self.O = ['公园散步', '购物', '清理房间']


    # 维特比算法
    def viterbi(self):
        # 初始化
        T = len(self.O)
        N = len(self.states)
        delta = np.zeros((T, N))
        fai = np.zeros((T, N), dtype=int)

        # 计算初始状态的delta
        delta[0] = self.pai * self.B[:, self.observations[self.O[0]]]

        # 递推计算delta和fai
        for t in range(1, T):
            for j in range(N):
                temp_delta = delta[t - 1] * self.A[:, j]
                max_delta_index = np.argmax(temp_delta)
                delta[t][j] = temp_delta[max_delta_index] * self.B[j][self.observations[self.O[t]]]
                fai[t][j] = max_delta_index

        # 回溯找到最优路径
        path = np.zeros(T, dtype=int)
        path[T - 1] = np.argmax(delta[T - 1])
        for t in range(T - 2, -1, -1):
            path[t] = fai[t + 1][path[t + 1]]
p
        # 返回最优路径和对应的天气状态
        weather_sequence = []
        for p in path:
            key = list(self.states.keys())[p]
            weather_sequence.append(key)
        return weather_sequence,delta


if __name__ == '__main__':
    HMM=HMM()
    # 调用维特比算法
    weather_sequence,delta = HMM.viterbi()

    # 输出结果
    print("前天天气:", weather_sequence[0])
    print("昨天天气:", weather_sequence[1])
    print("今天天气:", weather_sequence[2])
    print("概率矩阵为:\n",delta)

3.3 answer

在这里插入图片描述
  用这个去guess你的女朋友吧~~~~~~~

  • 3
    点赞
  • 0
    收藏
    觉得还不错? 一键收藏
  • 打赏
    打赏
  • 0
    评论
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包

打赏作者

幻兒

你的鼓励将是我创作的最大动力

¥1 ¥2 ¥4 ¥6 ¥10 ¥20
扫码支付:¥1
获取中
扫码支付

您的余额不足,请更换扫码支付或充值

打赏作者

实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值