{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "### 一.基本思路：马尔可夫链收敛到平稳态\n", "在[《12_05_PGM_马尔科夫链_初探及代码实现》](https://nbviewer.jupyter.org/github/zhulei227/ML_Notes/blob/master/notebooks/12_05_PGM_%E9%A9%AC%E5%B0%94%E7%A7%91%E5%A4%AB%E9%93%BE_%E5%88%9D%E6%8E%A2%E5%8F%8A%E4%BB%A3%E7%A0%81%E5%AE%9E%E7%8E%B0.ipynb)这一节，我们首次探索了马尔可夫链，并在本末讨论并说明了它的一个重要性质：**马尔可夫链的平稳态**，而且在后面的PageRank算法中使用了它的这个性质，这部分内容其实是要求这样一个问题： \n", "\n", "
**已知马尔可夫模型的状态转移概率矩阵$P$，求它的平稳态$\\pi^*$，使得$P\\pi^*=\\pi^*$**
\n", "\n", "将这个问题反过来，那就是**MCMC**（马尔科夫蒙特卡洛抽样法）的主要思想咯： \n", "\n", "
**已知某分布$\\pi$，求一个马尔可夫模型的状态转移概率矩阵$P^*$，使得$P^*\\pi=\\pi$**
\n", "\n", "于是，我们求一个$P^*$满足上面的条件即可，但是要用于MCMC，对$P^*$的要求更加苛刻一些：要保证在$P^*$的条件下，**平稳态必须是唯一**的，换言之，给定一个概率转移矩阵，它的平稳态可能有多个，比如下面的例子" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import pandas as pd\n", "import os\n", "os.chdir('../')\n", "from ml_models.pgm import SimpleMarkovModel" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "smm=SimpleMarkovModel(status_num=3)\n", "smm.P=np.asarray([\n", " [1,1/3,0],\n", " [0,1/3,0],\n", " [0,1/3,1]\n", "])" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
5010050010002000
00.4333330.5499945.500000e-015.500000e-015.500000e-01
10.2333330.0000129.750689e-251.358228e-482.635403e-96
20.3333330.4499944.500000e-014.500000e-014.500000e-01
\n", "
" ], "text/plain": [ " 50 100 500 1000 2000\n", "0 0.433333 0.549994 5.500000e-01 5.500000e-01 5.500000e-01\n", "1 0.233333 0.000012 9.750689e-25 1.358228e-48 2.635403e-96\n", "2 0.333333 0.449994 4.500000e-01 4.500000e-01 4.500000e-01" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "#定义初始状态\n", "init_prob=np.asarray([[0.2],[0.7],[0.1]])\n", "pd.DataFrame({50:smm.predict_prob_distribution(1,set_init_prob=init_prob).reshape(-1).tolist(),\n", " 100:smm.predict_prob_distribution(10,set_init_prob=init_prob).reshape(-1).tolist(),\n", " 500:smm.predict_prob_distribution(50,set_init_prob=init_prob).reshape(-1).tolist(),\n", " 1000:smm.predict_prob_distribution(100,set_init_prob=init_prob).reshape(-1).tolist(),\n", " 2000:smm.predict_prob_distribution(200,set_init_prob=init_prob).reshape(-1).tolist()})" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
5010050010002000
00.8333330.8499998.500000e-018.500000e-018.500000e-01
10.0333330.0000021.392956e-251.940325e-493.764862e-97
20.1333330.1499991.500000e-011.500000e-011.500000e-01
\n", "
" ], "text/plain": [ " 50 100 500 1000 2000\n", "0 0.833333 0.849999 8.500000e-01 8.500000e-01 8.500000e-01\n", "1 0.033333 0.000002 1.392956e-25 1.940325e-49 3.764862e-97\n", "2 0.133333 0.149999 1.500000e-01 1.500000e-01 1.500000e-01" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "#定义另一组初始状态\n", "init_prob=np.asarray([[0.8],[0.1],[0.1]])\n", "pd.DataFrame({50:smm.predict_prob_distribution(1,set_init_prob=init_prob).reshape(-1).tolist(),\n", " 100:smm.predict_prob_distribution(10,set_init_prob=init_prob).reshape(-1).tolist(),\n", " 500:smm.predict_prob_distribution(50,set_init_prob=init_prob).reshape(-1).tolist(),\n", " 1000:smm.predict_prob_distribution(100,set_init_prob=init_prob).reshape(-1).tolist(),\n", " 2000:smm.predict_prob_distribution(200,set_init_prob=init_prob).reshape(-1).tolist()})" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
5010050010002000
00.5666670.6499966.500000e-016.500000e-016.500000e-01
10.1666670.0000086.964778e-259.701626e-491.882431e-96
20.2666670.3499963.500000e-013.500000e-013.500000e-01
\n", "
" ], "text/plain": [ " 50 100 500 1000 2000\n", "0 0.566667 0.649996 6.500000e-01 6.500000e-01 6.500000e-01\n", "1 0.166667 0.000008 6.964778e-25 9.701626e-49 1.882431e-96\n", "2 0.266667 0.349996 3.500000e-01 3.500000e-01 3.500000e-01" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "#再定义另一组初始状态\n", "init_prob=np.asarray([[0.4],[0.5],[0.1]])\n", "pd.DataFrame({50:smm.predict_prob_distribution(1,set_init_prob=init_prob).reshape(-1).tolist(),\n", " 100:smm.predict_prob_distribution(10,set_init_prob=init_prob).reshape(-1).tolist(),\n", " 500:smm.predict_prob_distribution(50,set_init_prob=init_prob).reshape(-1).tolist(),\n", " 1000:smm.predict_prob_distribution(100,set_init_prob=init_prob).reshape(-1).tolist(),\n", " 2000:smm.predict_prob_distribution(200,set_init_prob=init_prob).reshape(-1).tolist()})" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "由上面的例子可以发现，由于初始状态的不同，平稳态可能有多个，所以我们构造的$P$必须要保证平稳态能唯一的收敛到我们的目标分布$\\pi$，不然，如上面的例子，假如我们的目标分布是$\\pi=[0.65,0,0.35]$，然后我们构造了一个概率转移矩阵$P=\\begin{bmatrix}\n", "1 & 1/3 & 0\\\\ \n", "0 & 1/3 &0 \\\\ \n", "0 & 1/3 & 1\n", "\\end{bmatrix}$ ，满足$P\\pi=\\pi$，然后，我们随便定义了一个初始点$x_0=[0.8,0.1,0.1]$，在$P$上随机游走采样$\\{x_1,x_2,...,x_n\\}$，最后发现它成功收敛到了$x_n\\rightarrow [0.85,0,0.15]$，哈哈哈哈哈~~~~ \n", "\n", "保证马尔科夫模型仅收敛到唯一平稳态的**充分条件**是有的！那就是**细致平衡方程**" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 二.保证平稳态唯一：细致平衡方程\n", "细致平衡方程的定义如下，设有马尔科夫链$X=\\{X_0,X_1,...,X_t,...\\}$，状态空间为$S$，转移概率矩阵为$P$，如果有状态分布$\\pi=(\\pi_1,\\pi_2,...)^T$，对任意时刻状态$i,j\\in S$，对任意时刻$t$均满足： \n", "\n", "$$\n", "P(X_t=i\\mid X_{t-1}=j)\\pi_j=P(X_{t-1}=j\\mid X_t=i)\\pi_i,i,j=1,2,....\n", "$$ \n", "\n", "或者简写为： \n", "\n", "$$\n", "p_{ij}\\pi_j=p_{ji}\\pi_i,i,j=1,2,....\n", "$$ \n", "\n", "为了更加直观，也可以写作： \n", "\n", "$$\n", "p(j\\rightarrow i)\\pi_j=p(i\\rightarrow j)\\pi_i\n", "$$\n", "这便是细致平衡方程，此时的马尔科夫链称为可逆马尔科夫链，显然对$P$矩阵按列求和为1，即$\\sum_{i}p_{ij}=1$或者$\\sum_ip(j\\rightarrow i)=1$ ，接下来简单证明一下，满足细致平衡条件，则有$P\\pi=\\pi$成立： \n", "\n", "$$\n", "(P\\pi)_i=\\sum_jp_{ij}\\pi_j=\\sum_jp_{ji}\\pi_i=\\pi_i\\sum_{j}p_{ji}=\\pi_i,i=1,2,...\n", "$$ \n", "\n", "**高维/连续状态空间**：上面只是对一维离散状态空间做的说明，而细致平衡方程对高维空间或者连续状态空间一样是成立的，只需对相应符号做修改即可，比如将$p_{ij}$修改为一个函数的形式$p(状态j,状态i)$，求和符号$\\sum$需要替换为积分符号$\\int$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 三.小结\n", "最后小结一下，这一节的重点引出细致平衡方程： \n", "\n", "$$\n", "p_{ij}\\pi_j=p_{ji}\\pi_i\n", "$$\n", "\n", "因为在满足该方程约束的条件下，我们构造的$P^*$，必然有$P^*\\pi=\\pi$，且$\\pi$唯一（这里$\\pi$是我们的目标分布），这也是检验MCMC是否有效的一条重要标准，下一节将要介绍的Metropolis-Hastings算法（MH算法）便是满足该约束条件的一套算法\n", "\n", "#### 采样步骤\n", "\n", "再对MCMC的采样步骤做个整理 \n", "\n", "（1）首先在随机变量$x$的状态空间上构造一个马尔科夫链$P$，使它能平稳且唯一的收敛到我们的目标分布$p(x)$； \n", "\n", "（2）从状态空间某一点$x_0$出发，用构造的马尔科夫链$P$进行随机游走，产生样本序列$x_0,x_1,...,x_t,....$； \n", "\n", "（3）确定一个正整数$m,n$（$m