YidaXu-ML/sampling-methods-part2.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 采样方法系列2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## MCMC\n",
    "\n",
    "这一节将介绍MCMC（Markov Chain Monte Carlo）采样算法，MCMC是一种用马氏链从随机分布取样的算法，之前步骤的作为底本。步数越多，结果越好。首先我将介绍MH算法，接着介绍Gibbs采样，然后证明Gibbs采样是MH算法的一种特例。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 马尔科夫链\n",
    "\n",
    "本节以一个简单的例子来回顾马氏过程中的一些基本概念。社会学家经常把人按其经济状况分成 3 类：下层 (lower-class)、中层 (middle-class)、上层 (upper-class)，我们用 1、2、3 分别代表这三个阶层。社会学家们发现决定一个人的收入阶层的最重要的因素就是其父母的收入阶层。如果一个人的收入属于下层类别，那么他的孩子属于下层收入的概率是 0.65，属于中层收入的概率是 0.28，属于上层收入的概率是 0.07。事实上，从父代到子代，收入阶层的变化的转移概率如下\n",
    "\n",
    "|-|\t下层\t|中层|上层|\n",
    "|:|:|:|:|\n",
    "|下层|0.65|0.28|0.07|\n",
    "|中层|0.15|0.67|0.18|\n",
    "|上层|0.12|0.36|0.52|\n",
    "\n",
    "它们之间的概率转移图为：\n",
    "![](https://uploads.cosx.org/2013/01/markov-transition.png)\n",
    "\n",
    "现在，我们要计算第N代后，整个社会中各阶层人数的比例：\n",
    "\n",
    "计算概率转移矩阵\n",
    "\n",
    "$$P =\\left[ \\begin{matrix}  0.65 &  0.28 & 0.07\\\\  0.15 &  0.67 & 0.18 \\\\  0.12 &  0.36 & 0.52 \\end{matrix} \\right]$$\n",
    "\n",
    "\n",
    "假设当前这一代人处在下层、中层、上层的人的比例是概率分布向量$\\pi_0=[\\pi_0(1), \\pi_0(2), \\pi_0(3)]$，那么他们的子女的分布比例将是$\\pi_1=\\pi_0P$，他们的孙子代的分布比例将是$\\pi_2 = \\pi_1P=\\pi_0P^2$，……, 第$n$代子孙的收入分布比例将是 $\\pi_n = \\pi_{n-1}P = \\pi_0P^n$\n",
    "\n",
    "接下来我们用程序来模拟这个过程，首先定义概率转移矩阵$P$和子孙代数$N$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from numpy import linalg as LA\n",
    "\n",
    "def compute_dist(init_dist, P, N):\n",
    "    '''\n",
    "    计算子代分布\n",
    "    init_dist: 初始分布\n",
    "    P: 概率转移矩阵\n",
    "    N: 子孙代数\n",
    "    '''\n",
    "    assert P.shape[0] == P.shape[1]\n",
    "    for n in range(N + 1):\n",
    "        current_dist = np.dot(init_dist, LA.matrix_power(P, n))\n",
    "        print '第{}代人'.format(n), current_dist \n",
    "    return current_dist\n",
    "\n",
    "# 概率转移矩阵\n",
    "P = np.array([[.65, .28, .07],[.15,.67,.18],[.12,.36,.52]])\n",
    "# 子孙代数\n",
    "N = 10\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "第一次实验，假设当代人处在下层、中层、上层的比例是$\\pi_0=[0.51,0.26,0.23]$，那么让我们来看一下$N=10$代以后，他们的子孙在三个阶层的分布分别是多少"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "第0代人 [ 0.51  0.26  0.23]\n",
      "第1代人 [ 0.3981  0.3998  0.2021]\n",
      "第2代人 [ 0.342987  0.45209   0.204923]\n",
      "第3代人 [ 0.31534581  0.47270894  0.21194525]\n",
      "第4代人 [ 0.30131455  0.48131211  0.21737335]\n",
      "第5代人 [ 0.29413607  0.48510159  0.22076234]\n",
      "第6代人 [ 0.29044517  0.48685061  0.22270423]\n",
      "第7代人 [ 0.28854146  0.48768807  0.22377047]\n",
      "第8代人 [ 0.28755761  0.48809999  0.2243424 ]\n",
      "第9代人 [ 0.28704854  0.48830639  0.22464508]\n",
      "第10代人 [ 0.28678492  0.4884111   0.22480399]\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "array([ 0.28678492,  0.4884111 ,  0.22480399])"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "init_state1 = np.array([.51, .26, .23])\n",
    "compute_dist(init_state1, P, N)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "第二次实验，我们试着将当代人在下层、中层、上层的比例固定为$\\pi_0=[0.32,0.47,0.21]$，看看会发生什么"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "第0代人 [ 0.32  0.47  0.21]\n",
      "第1代人 [ 0.3037  0.4801  0.2162]\n",
      "第2代人 [ 0.295364  0.484535  0.220101]\n",
      "第3代人 [ 0.29107897  0.48657673  0.2223443 ]\n",
      "第4代人 [ 0.28886916  0.48755247  0.22357838]\n",
      "第5代人 [ 0.28772723  0.48803173  0.22424104]\n",
      "第6代人 [ 0.28713638  0.48827166  0.22459196]\n",
      "第7代人 [ 0.28683043  0.4883933   0.22477626]\n",
      "第8代人 [ 0.28667193  0.48845549  0.22487258]\n",
      "第9代人 [ 0.28658979  0.48848745  0.22492277]\n",
      "第10代人 [ 0.28654721  0.48850393  0.22494886]\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "array([ 0.28654721,  0.48850393,  0.22494886])"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "init_state2 = np.array([.32, .47, .21])\n",
    "compute_dist(init_state2, P, N)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "我们发现虽然两次实验的初始分布不同，但最终经过若干代后，三个阶层的比例都会最终固化。实际上，只要$N$足够大，无论怎么初始化起始分布，最终都会收敛于某个分布，我们把这个分布叫做平稳分布。这个性质不是该马氏链独有的，绝大部分马氏链都有这个性质。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 细致平衡条件\n",
    "\n",
    "\n",
    "细致平衡条件是MCMC算法能够工作的前提。  \n",
    "\n",
    ">细致平衡条件  对于定义在状态空间$\\mathcal{X}$上的概率密度函数$\\pi$，如果对于任意的$i,j\\in \\mathcal{X}$，满足$\\pi_i P_{ij}=\\pi_j P_{ji}$，那么称$\\pi$关于$P$满足细致平衡条件（Detailed Balance）。\n",
    "\n",
    "其中，$P$是马氏链的转移概率矩阵，$P_{ij}$表示从状态$i$跳转到状态$j$的概率$P_{ij}=P(x^{(t)}=j|x^{(t-1)}=i)$。  \n",
    "从细致平衡条件可以推导出许多有用的结论  \n",
    "\n",
    "* 结论一  \n",
    "$\\pi$是$P$对应马氏链的平稳分布  \n",
    "证明：  \n",
    "\n",
    "$$ \\sum_i \\pi_i P_{ij}=\\sum_i \\pi_j P_{ji}=\\pi_j \\sum_i P_{ji}=\\pi_j$$\n",
    "\n",
    "由此可知，如果一个分布满足细致平衡条件，那么它一定是某个马氏链的平稳分布。我们可以通过构造马氏链的转移概率矩阵$P$，使得该马氏链的平稳分布刚好是要采样的分布。\n",
    "\n",
    "* 结论二  \n",
    "满足细致平衡条件的马氏链称为可逆的马氏链（reversible Markov chain) ，其性质如下：  \n",
    "如果$\\pi,P$满足细致平衡，且序列$(x_0,...,x_n)$服从以$\\pi,P$为参数的马氏链 $MC(\\pi, P)$，那么其逆序列$(x_n,...,x_0)$也服从该马氏链\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Metropolis–Hastings算法\n",
    "\n",
    "Metropolis–Hastings算法通过构造一个马氏链，并使其平稳分布为被采样分布，通过在该马氏链的平稳分布采样就获得了目标分布的样本。\n",
    "\n",
    "Metropolis–Hastings算法将细致平衡条件中的转移概率$P(x^*|x)$分解为提议概率$q(x^*|x)$和接受概率$A(x^*|x)$：\n",
    "$$P(x^*|x) = q(x^*|x) A(x^*|x)$$\n",
    "\n",
    "提议概率定义为$x$跳转到$x^*$的概率；接受概率表示接受$x^*$的概率。\n",
    "\n",
    "于是\n",
    "\n",
    "$$A(x^*|x)=\\frac{P(x^*|x)}{q(x^*|x)}$$\n",
    "\n",
    "定义接受率为\n",
    "\n",
    "$$ \\alpha(x^*|x) = \\frac{A(x^*|x)}{A(x|x^*)}=\\frac{P(x^*|x)q(x|x^*)}{P(x|x^*)q(x^*|x)}=\\frac{\\pi(x^*)q(x|x^*)}{\\pi(x)q(x^*|x)}$$\n",
    "\n",
    "这里用到了细致平衡条件\n",
    "\n",
    "$$\\pi(x)P(x^*|x)=\\pi(x^*)P(x|x^*)\\rightarrow \\frac{P(x^*|x)}{P(x|x^*)}=\\frac{\\pi(x^*)}{\\pi(x)}$$\n",
    "\n",
    "接受率的含义是从$x$转移到$x^*$的接受概率与从$x^*$转移到$x$的接受概率的比值，如果接受率大于1，则表明样本点更有可能出现在$x^*$，选择接受；如果接受率不大于1，则表明有$\\alpha$的可能性接受，$1-\\alpha$的可能性拒绝。最终我们倾向于留在高概率密度的地方，然后仅偶尔跑到低概率状态的地方（这也就是MH算法直观上的运行机理）。 \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "\n",
    "Metropolis–Hastings算法的流程如下图所示：\n",
    "![](http://7xikew.com1.z0.glb.clouddn.com/Metropolis-Hasting-Algorithm.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "其中$\\alpha(x^*|x)=\\min \\big(1, \\frac{\\pi(x^*)q(x|x^*)}{\\pi(x)q(x^*|x)}\\big)$是接受率（accept rate），主题到MH算法给接受率加了一个上界。\n",
    "\n",
    "简单描述一下算法的工作流程：  \n",
    "1）首先初始化马氏链的其实状态$x^{(0)}$  \n",
    "2）从$(0,1)$之间的均匀分布采样$u$，从转移概率分布$q(x^*|x^{(i)})$采样$x^*$  \n",
    "3）如果$u$小于采样率$\\alpha(x^*|x)$，则接受$x^*$；否则拒绝$x^*$  \n",
    "\n",
    "## MH算法细致平衡条件证明\n",
    "MH算法之所以能work，在于它满足了细致平衡条件。\n",
    "\n",
    "证明：  \n",
    "\n",
    "根据上一节推导的的结论一，我们只要证明$\\pi(x)$满足细致平衡条件。\n",
    "\n",
    "将马氏链的转移概率分解为\n",
    "$$P(x^*|x)=q(x^*|x)A(x^*|x)$$  \n",
    "\n",
    "$P(x^*|x)$由两个步骤构成：\n",
    "\n",
    "1）从提议分布$q(x^*|x)$产生$x^*$  \n",
    "2）以概率$A(x^*|x)$接受$x^*$\n",
    "\n",
    "\n",
    "接下来证明细致平衡条件\n",
    "\n",
    "$$\\begin{aligned}\\pi(x)P(x^*|x)&=\\pi(x)q(x^*|x)A(x^*|x)\\\\&=\\pi(x)q(x^*|x)\\min \\big(1, \\frac{\\pi(x^*)q(x|x^*)}{\\pi(x)q(x^*|x)}\\big)\\\\&=\\min(\\pi(x)q(x^*|x),\\pi(x^*)q(x|x^*))\\\\&=\\pi(x^*)q(x|x^*)\\min \\big(1, \\frac{\\pi(x)q(x^*|x)}{\\pi(x^*)q(x|x^*)}\\big)\\\\&=\\pi(x^*)q(x|x^*) A(x|x^*)\\\\&=\\pi(x^*)P(x|x^*)\\end{aligned}$$\n",
    "\n",
    "\n",
    "由此，我们可知我们只要构造一个转移概率分布为$P(x^*|x)$的马氏链，并在该马氏链的平稳分布上进行抽样，就可以获得到目标分布$\\pi(x)$的样本。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Gibbs采样算法\n",
    "\n",
    "给定一个起始样本点$(x^{(1)},y^{(1)},z^{(1)})^\\top$，我们想从目标分布$p(x,y,z)$采样得到$\\{(x^{(2)},y^{(2)},z^{(2)})^\\top,(x^{(3)},y^{(3)},z^{(3)})^\\top,...,(x^{(n)},y^{(n)},z^{(n)})^\\top\\}$，Gibbs采样算法的工作流程是这样的：  \n",
    "采样$(x^{(2)},y^{(2)},z^{(2)})^\\top$的过程：\n",
    "$$ x^{(2)}\\sim p(x|y^{(1)},z^{(1)})\\\\ y^{(2)}\\sim p(y|x^{(2)},z^{(1)})\\\\ z^{(2)}\\sim p(z|x^{(2)},y^{(2)})$$\n",
    "采样$(x^{(3)},y^{(3)},z^{(3)})^\\top$的过程： \n",
    "$$ x^{(3)}\\sim p(x|y^{(2)},z^{(2)})\\\\ y^{(3)}\\sim p(y|x^{(3)},z^{(2)})\\\\ z^{(3)}\\sim p(z|x^{(3)},y^{(3)})$$\n",
    "\n",
    "这个过程理解起来很简单，循环遍历$x,y,z$三个坐标轴，每次固定其他两个变量（比如$y,z$），然后对剩余的一个变量（比如$x$）进行采样。  \n",
    "\n",
    "考虑更一般的情况，假设要采样的随机向量有$K$维，定义$x_{-i}=\\{x_1,...,x_{i-1},x_{i+1},...,x_K\\}$，表示去掉$x_i$后其他维度上的随机变量构成的集合。\n",
    "\n",
    "更一般情况下，Gibbs采样算法定义如下：  \n",
    "> n维Gibbs采样算法  \n",
    "1.initilize $(x^{(0)}_1,...,x^{(0)}_K)$  \n",
    "2.for t in 0,1,2,...  \n",
    "&emsp;&emsp;$x_1^{(t+1)}\\sim p(x_1|x_2^{(t)},...,x_K^{(t)})$  \n",
    "&emsp;&emsp;$x_2^{(t+1)}\\sim p(x_2|x_1^{(t+1)},x_3^{(t)},...,x_K^{(t)})$  \n",
    "&emsp;&emsp;...  \n",
    "&emsp;&emsp;$x_i^{(t+1)}\\sim p(x_i|x_1^{(t+1)},...,x_{i-1}^{(t+1)},x_{i+1}^{(t)},...,x_K^{(t)})$  \n",
    "&emsp;&emsp;$x_K^{(t+1)}\\sim p(x_K|x_1^{(t+1)},...,x_{K-1}^{(t+1)})$\n",
    "\n",
    "其中\n",
    "$$p(x_i|x_{-i})=p(x_i|x_1,...,x_{i-1},x_{i+1},...,x_K)=\\frac{p(x_1,...,x_K)}{p(x_1,...,x_{i-1},x_{i+1},...,x_K)}$$\n",
    "是Gibbs采样的提议分布，也称为完全条件概率"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Gibbs采样细致平衡条件证明  \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "这个问题可以转化为证明Gibbs采样是MH算法的一种特殊情况。假设我们要对$x_i$进行采样，考察接受率\n",
    "$$\\alpha(x^*|x)=\\frac{\\pi(x^*)q(x|x^*)}{\\pi(x)q(x^*|x)}=\\frac{p(x^*)p(x_i|x_{-i}^*)}{p(x)p(x_i^*|x_{-i})}=\\frac{p(x_i^*|x^*_{-i})p(x^*_{-i})p(x_i|x_{-i}^*)}{p(x_i|x^{-i})p(x_{-i})p(x_i^*|x_{-i})}$$\n",
    "\n",
    "注意到$x^*_{-i}=x_{-i}$，于是\n",
    "\n",
    "$$\\alpha(x^*|x)=\\frac{p(x^*|x^*_{-i})p(x^*_{-i})p(x_i|x_{-i}^*)}{p(x|x^{-i})p(x_i)p(x_i^*|x_{-i})}=\\frac{p(x_i^*|x_{-i})p(x_{-i})p(x_i|x_{-i})}{p(x_i|x_{-i})p(x_{-i})p(x_i^*|x_{-i})}=1$$\n",
    "\n",
    "由此我们知道，Gibbs采样是MH算法在接受率等于1时的特殊情形，也就是说Gibbs采样的过程只要一直接受提议分布产生的样本点即可。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Gibbs采样的并行化\n",
    "\n",
    "未完待续"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "参考：  \n",
    "1.[LDA-math-MCMC 和 Gibbs Sampling by 靳志辉](https://cosx.org/2013/01/lda-math-mcmc-and-gibbs-sampling)  \n",
    "2.[LDA算法理解](https://www.zybuluo.com/Dounm/note/435982)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}