《数据挖掘》第六课笔记

---

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

  在计算积分值时，蒙特拉洛的主要内容就是从假设分布中抽取n次，估计积分值。对于概率分布不是常见的分布，一个可行的办法是采用接受-拒绝采样来得到该分布的样本。既然 $p(x)$太复杂在程序中没法直接采样，那么我设定一个程序可采样的分布$q(x)$ 比如正态分布，然后按照一定的方法拒绝某些样本，以达到接近 $p(x)$ 分布的目的。

  具体采用过程如下，设定一个方便采样的常用概率分布函数 $q(x)$，以及一个常量 $k$，使得 $p(x)$总在 $kq(x)$ 的下方，如上图。
  首先，采样得到$q(x)$的一个样本$z_0$，然后，从均匀分布$(0,kq(z_0))$中采样得到一个值$u$。如果$u$落在了上图中的灰色区域，则拒绝这次抽样，否则接受这个样本$z_0$。重复以上过程得到$n$个接受的样本 $z_0,z_1,\dots z_{n-1}$，则最后的蒙特卡罗求解结果为：$$\frac{p(z_i)}{kq(z_i)}$$  整个过程中，我们通过一系列的接受拒绝决策来达到用$q(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$$\ge$$0$，及状态$i,j,i_0,i_1,L,i_{n-1}$，有（马氏性）$$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)$$  即马尔可夫链 { X n , n ≥ 0 } \{ {X_{\rm{n}}},n \ge 0\} {Xn​,n≥0}的有限维分布完全由初始分布 P { X 0 = i 0 } P\{ {X_0} = {i_0}\} P{X0​=i0​}和条件概率（一步转移概率）$P(X_n=j|X_{n-1}=i)$确定

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

2.1 状态的分类

  定义1:若存在$n$$\ge$$0$使得$P_{ij}^{(n)}$$>$$0$,则称从状态$i$可达状态$j$，记$i\to j$。如果$i\to j$且$j\to i$,则称$i$和$j$互通，记$i\leftrightarrow j$
  将任何两个互通的状态归为一类：（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:n≥1,Pii(n)​>0}=∅，则称它的最大公约数$d=d(j)$为状态$i$的周期，若 { d = 1 , 称 i 为非周期的 d > 1 , 称 i 为周期的 \color{red}\{ _{d = 1,称i为非周期的}^{d > 1,称i为周期的} {d=1,称i为非周期的d>1,称i为周期的​特.若上述集合为空集，则称$i$的周期为无穷大
  定理:若状态$i,j$同属一类，则$d(i)=d(j)$

2.3 平稳分布（不变分布）与极限分布

  平稳分布:设马氏链 { X n , n ≥ 0 } \{ {X_{\rm{n}}},n \ge 0\} {Xn​,n≥0}有转移概率矩阵$P=({\mathrm{P}}_{\mathrm{ij}})$，若存在一个概率分布 { π i , i ∈ S } \{ {\pi _i},i \in S\} {πi​,i∈S}，满足${\pi }_j=\sum _{i\in S}{\pi }_i{P}_{ij}$，则称 { π i , i ∈ S } \color{red}\{ {\pi _i},i \in S\} {πi​,i∈S}为该链的平稳分布（不变分布）
  分析：
   令$\pi =({\pi }_1,{\pi }_2,\dots )$，则上式为$\pi =\pi P$
   (1).由$\pi =\pi P$可知，$\pi (I-P)=0$，故$I$是矩阵$P$的左特征值，平稳分布$\pi$是$P$的左特征向量
   (2).两边同乘$P$，得$\pi =\pi P=\pi P^2=\dots =\pi P^n$

  极限分布:称Markov链是遍历的，如果所有状态相通且均是周期为1的正常返状态，对于遍历的Markov链，极限$\underset{n\to \infty }{\lim }{P}_{ij}^{(n)}={\pi }_j=\frac{1}{u_j}$，$j\in S$称为Markov链的极限分布

3.基于马尔可夫链采样

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

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

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

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

---

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

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

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

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

（1）概率计算问题
  给定：$\lambda =(A,B,\pi )$  $O=(o_1,o_2,\cdots ,o_T)$
  计算：$P(O|\lambda )$
（2）学习问题
  已知：$O=(o_1,o_2,\cdots ,o_T)$
  估计：$\lambda =(A,B,\pi )$，使$P(O|\lambda )$最大
（3）预测问题(解码)
  已知：$\lambda =(A,B,\pi )$  $O=(o_1,o_2,\cdots ,o_T)$
  求:使$P(O|\lambda )$最大的状态序列$I=(i_1,i_2,\cdots ,i_T)$

3.隐马氏模型的应用

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

---

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

1.暴力求解

  1.一个简单问题
  知道骰子有几种，每种骰子是什么，每次掷的都是什么骰子，根据掷骰子掷出的结果，求产生这个结果的概率。

$P=P(D6)\ast P(D6\to 1)\ast P(D6\to D8)\ast P(D8\to 6)\ast P(D8\to D8)\ast P(D8\to 3)$
$=\frac{1}{3}\ast \frac{1}{6}\ast \frac{1}{3}\ast \frac{1}{8}\ast \frac{1}{3}\ast \frac{1}{8}$

  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)$

2.前向算法

  仅做一个基本介绍，前向算法的时间复杂度是$O(N^{2}T)$
  前向概率定义：给定隐马尔可夫模型
定义到时刻$t$部分观测序列：$o_1,o_2,\cdots ,o_i$，且状态为$q_i$的概率为前向概率，记作${\alpha }_t(i)=P(o_1,o_2,\cdots ,o_t,i_t=q_i|\lambda )$
  算法（观测序列概率的前向算法）
  输入：隐马尔可夫模型$\lambda$，观测序列$O$
  输出：观测序列概率$P(O|\lambda )$
  初值：${\alpha }_t(i)={\pi }_i{b}_{\mathrm{i}}(o_1)$，$i=1,2,...,N$
  递推：${\alpha }_{t+1}(i)=\left[\sum _{i=1}^N{\alpha }_t(j){\alpha }_{ji}\right]b_i(o_{t+1})$，$i=1,2,...,N$
  终止：$P(O|\lambda )=\sum _{i=1}^N{\alpha }_T(i)$

---

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

1.算法（维特比算法）

  输入：模型$\lambda =(A,B,\pi )$和观测$O=(o_1,o_2,\cdots ,o_T)$
  输出：最优路径$I^{\ast }=(i_1^{\ast },i_2^{\ast },\cdots ,i_T^{\ast })$
  （1）初始化$${\delta }_1(i)={\pi }_i{b}_i(o_1)，i=1,2,...N$$$${\varphi }_1(i)=0，i=1,2,...N$$  （2）递推.对$t=2,3,...,T$$${\delta }_t(i)=\underset{1\le j\le N}{\max }[{\delta }_{t-1}(j)a_{ji}]b_i(o_t)，i=1,2,...N$$$${\varphi }_t(i)=\arg \underset{1\le j\le N}{\max }[{\delta }_{t-1}(j)a_{ji}]，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你的女朋友吧~~~~~~~