马尔科夫链蒙特卡洛应用初窥

在本文中，笔者会持续更新一些应用实例来帮助理解MCMC，以及一些现有的工具包使得我们更快的上手。(我的这句中文为嘛看起来像机翻) 用采样的方法不单单可以帮我们找到比较好的参数，也可以帮我们发现这些参数的分布。

常用的工具

一个简单的例子

根据贝叶斯公式: $p(\theta | x) = \frac{p(\theta)p(x|\theta)}{p(x)} \propto p(\theta)p(x|\theta)$

其中$x$表示我们的观测数据, $p(\theta)$ 表示参数的先验分布，可以理解为在看到数据$x$之前参数的分布，$p(\theta | x)$表示的是后验分布，可以理解为在看到数据之后，参数的分布，$p(x |\theta)$是’likelihood’，表示在当前参数下，我们观测数据的概率。

“A Bayesian is one who, vaguely expecting a horse and catching a glimpse of a donkey, strongly concludes he has seen a mule” (Senn)

例子[1]: 体重和身高的表达式(一条样本):

\[weight_i = \beta_0 + \beta_1 height_i + e_i \quad e_i \sim N(0, \sigma^2)\]

其中各个参数的分布如下:

\[\beta_0 \sim N(0, m_0), \beta_1 \sim N(0, m_1), \sigma^2 \sim \Gamma^{-1}(\epsilon, \epsilon) \quad m_0=m_1=10^6, \epsilon=10^{-3}\]

我们的采样目标概率密度分布是后验分布: $p(\beta_0, \beta_1, \sigma^2 | y)$，但是如果直接在这个分布上采样会有多维的问题，所以我们可以用吉布斯采样，把它变成条件后验分布，分别采样:

\[p(\beta_0|y,\beta_1, \sigma^2),p(\beta_1|y, \beta_0, \sigma^2), p(\sigma^2|y, \beta_0, \beta_1)\]

吉布斯采样: 先初始化这几个参数，然后按照上面的条件后验分布分别为各个参数采样，将三个新参数加入到样本集，继续这个循环进行采样

在之前的吉布斯采样中，我们知道因为转移函数也是目标概率密度函数$p$，所以我们希望我们的$p$能比较容易抽样。所以我们来分别看看这几个后验分布表达式(省略计算过程，详细请见[1]):

\[\begin{align}p(\beta_0|y,\beta_1, \sigma^2) &\propto p(\beta_0)p(y|\beta_0, \beta_1, \sigma^2)\\ m_0 &\rightarrow \infty, \beta_0 \sim N(\frac{1}{N}\sum_i (y_i-x_i\beta_1), \frac{\sigma^2}{N})\end{align}\] \[\begin{align}p(\beta_1|y,\beta_0, \sigma^2) &\propto p(\beta_1)p(y|\beta_0, \beta_1, \sigma^2)\\ m_0 &\rightarrow \infty, \beta_1 \sim N(\frac{\sum_i x_iy_i-\beta_0\sum_ix_i}{\sum_ix_i^2}, \frac{\sigma^2}{\sum_ix_i^2})\end{align}\] \[\begin{align}p(1/\sigma^2|y,\beta_0, \beta_1) &\propto p(1/\sigma^2)p(y|\beta_0, \beta_1, \sigma^2)\\ p(1/\sigma^2|y,\beta_0, \beta_1) &\sim \Gamma(a,b) \\ &where \quad a=\epsilon +\frac{N}{2}, b=\epsilon + \frac{1}{2}\sum_ie_i^2 \end{align}\]

我们可以发现这三个条件后验分布都是比较容易采样的函数。因此只要带入到我们之前的吉布斯采样函数中即可。

LDA

我们在这一节中来简单地看一下吉布斯采样如何求解潜在狄利克雷分配(latent Dirichlet allocation, LDA)，LDA的另一种求解方法是用变分EM算法，下篇文章我们来一探变分EM算法。下面我们先回顾几个知识点，一些证明会略掉，大家可以去李航老师的《统计学习方法》[2]中得到具体的证明过程。

多项式分布 若多元离散随机变量$X = (X_1, X_2, …, X_k)$的概率质量函数为：

\[\begin{align} P(X_1=n_1, X_2=n_2, \dots, Xk=n_k) &= \frac{n!}{n_1!n_2!\dots n_k!}p_1^{n_1}p_2^{n_2}\dots p_k{n_k}\\ &=\frac{n!}{\prod_{i=1}^k n_i!} \prod_{i=1}^kp_i^{n_i} \end{align}\]

其中 $ p=(p_1,p2,\dots, p_k), \quad p_i \ge 0, i=1,2, \dots, k, \quad \sum_{i=1}^k p_i=1, \sum_{i=1}^k n_i = n, $, 则称随机变量 $X$ 服从参数为 $(n,p)$ 的多项式分布，记作 $X \sim Mult(n,p)$ 。当实验的次数 $n=1$ ，多项式分布变为类别分布。

狄利克雷分布 （Dirichlet distribution），在贝叶斯学习中，狄利克雷分布常常作为多项式分布的先验分布使用(狄利克雷分布是多项式分布的共轭先验)。

\[\begin{align} p(\theta|\alpha) &= \frac{1}{B(\alpha)} \prod_{i=1}^k \theta_i^{\alpha_i-1} \\ B(\alpha) &= \frac{\prod_{i=1}^k \Gamma(\alpha_i)}{\Gamma(\sum_{i=1}^k \alpha_i)}\\ B(\alpha)&=\int\prod_{i=1}^k \theta_i^{\alpha_i-1}d\theta 是规范化因子，称为多元贝塔函数（证略） \end{align}\]

其中 $\sum_{i=1}^k \theta_i=1, \theta_i \ge0, \alpha=(\alpha_1, \alpha_2, \dots, \alpha_k), \alpha_i>0, i=1,2, \dots, k$ , 则称随机变量$\theta$服从参数为$\alpha$的狄利克雷分布，记作 $\theta \sim Dir(\alpha)$。$\Gamma(s)$为伽马函数。

接下来的这个性质在计算LDA的时候会经常用到所以也在这里贴一下。

设 $\mathcal{W} = {w1,w2,\dots,wk}$ 是由 $k$ 个元素组成的集合。随机变量 $X$ 服从$\mathcal{W}$ 上的多项式分布，$X \sim Mult(n,\theta)$ 其中 $n=(n_1, n_2, \dots, n_k)$ 和 $\theta=(\theta_1, \theta_2, \dots, \theta_k)$ 是参数。参数$n$为从 $\mathcal{W}$ 中重复独立抽取样本的次数，$n_i$ 为样本中 $w_i$ 出现的次数 $(i=1,2,\dots,k)$ ; 参数 $\theta_i$ 为 $w_i$ 出现的概率 $(i=1,2, \dots, k)$

将样本数据表示为 $D$，目标是计算在样本数据 $D$ 给定条件下参数 $\theta$ 的后验概率 $p(\theta | D)$。对于给定的样本数据 $D$，似然函数是:

\[p(D|\theta) = \theta_1^{n_1} \theta_2^{n_2} \dots \theta_k^{n_k} = \prod_{i=1}^k \theta_i^{n_i}\]

随机变量 $\theta$ 服从如前所述的狄利克雷分布, 记作: $p(\theta |\alpha) = Dir(\theta | \alpha), \alpha_i>0$

那么我们就可以得到我们的后验分布(具体计算省略):

\[\begin{align}p(\theta|D,\alpha) &= \frac{p(D|\theta)p(\theta|\alpha)}{p(D|\alpha)} \\ & \dots \\ &= Dir(\theta | \alpha+n) \end{align}\]

上述的性质为狄利克雷的后验分布的参数等于狄利克雷先验分布的参数 $\alpha=(\alpha_1, \alpha_2, \dots, \alpha_k)$ 加上多项式分布的观测计数 $n=(n_1, n_2, \dots, n_k)$, 因为就好像是实验之前就已经观察到了计数 $\alpha$，所以它也被称为先验伪计数。

LDA模型

终于到我们的LDA部分啦。我们首先来看一下我们期望解决的问题。通常我们有一堆文本，每个文本中有话题出现的概率不太一样，比如一篇关于体育的文章，比赛和羽毛球是比较可能的话题；每个话题下出现的单词概率分布也不同，比如话题比赛下，得分、苦战这样的词汇会更集中。这里，我们假设每一篇文本和所有的话题、每一个话题和所有的单词都是多项式分布，并且话题和单词都是狄利克雷先验分布。

这里我们定义我们的参数如下: 单词集合 $W= {w_1, \dots, w_v, \dots, w_V}$其中 $w_v$ 是第 $v$ 个单词, $V$ 是单词的总个数。文本集合$D={\boldsymbol{w_1, \dots, w_m, \dots, w_M }}$, 其中$w_m$是第m个文本， $M$是文本的个数。文本 $\boldsymbol{w_m}$ 是一个单词序列 $\boldsymbol{w_m}=(w_{m1}, \dots, w_{mn}, w_{mN_m})$ 其中 $w_{mn}$ 是文本 $\boldsymbol{w_m}$ 的第 $n$ 个单词，$N_m$ 是文本 $\boldsymbol{w_m}$ 中的单词个数。话题集合 $Z = { z_1, \dots, z_k, \dots, z_K}$, 其中 $z_k$ 是第 $k$ 个话题，$k=1,2, \dots, K$，$K$ 是话题个数。

每一个话题 $z_k$ 下的单词概率分布 $p(w |z_k)$ 服从多项式分布(类别分布)，其参数为 $\varphi_k$，它服从狄利克雷分布 (先验分布)，超参数为 $\beta = (\beta_1, \beta_2, \dots, \beta_V)$ 。参数 $ \varphi_k = (\varphi_{k1}, \varphi_{k2}, \dots, \varphi_{kV})$ 中的 $\varphi_k$ 表示话题 $z_k$ 生成单词 $w_v$ 的概率。$\boldsymbol{\varphi}={ \varphi_k}_{k=1}^K$ 组成了所有话题的参数矩阵:

\[\begin{pmatrix} \varphi_{11} & \varphi_{12} & \dots & \varphi_{1V} \\ \varphi_{21} & \varphi_{22} & \dots & \varphi_{2V} \\ \dots & \dots &\dots & \dots \\ \varphi_{k1} & \varphi_{k2}& \dots & \varphi_{kV} \end{pmatrix}\]

矩阵中的每一行向量都是从 $Dir(\beta)$ 中生成的。

同理，每一个文本 $\boldsymbol{w_m}$ 下的话题概率分布 $p(z |\boldsymbol{w_m})$ 服从多项式分布(类别分布)，其参数为 $\theta_m$，它服从狄利克雷分布 (先验分布)，超参数为 $\alpha= (\alpha_1, \alpha_2, \dots, \alpha_K)$ 。参数 $ \theta_m = (\theta_{m1}, \theta_{m2}, \dots, \theta_{mK})$ 中的 $\theta_{mk}$ 表示文本 $\boldsymbol{w_m}$ 生成话题 $z_k$ 的概率。$\boldsymbol{\theta}={ \theta_m}_{m=1}^M$ 组成了所有文本的参数矩阵:

\[\begin{pmatrix} \theta_{11} & \theta_{12} & \dots & \theta_{1K} \\ \theta_{21} & \theta_{22} & \dots & \theta_{2K} \\ \dots & \dots &\dots & \dots \\ \theta_{m1} & \theta_{m2}& \dots & \theta_{mK} \end{pmatrix}\]

矩阵中的每一行向量都是从 $Dir(\alpha)$ 中生成的。

首先我们分别用 $Dir(\beta)$ 和 $Dir(\alpha)$ 生成好每个话题的单词分布，和每个文本下的话题分布。接着对于某个文本$\boldsymbol{w_m}$，我们用其相应的参数 $\theta_m$ 按照多项式分布 $Multi(\theta_m)$ 随机生成一个话题 $z_{mn}$。找到这个话题 $z_{mn}$ 对应的单词分布参数 $\varphi_{z_{mn}}$，按照多项式分布 $Multi(\varphi_{z_{mn}})$ 随机生成一个单词 $w_{mn}$。重复这个过程，生成所有的文本单词。

LDA 参数学习

从上面的论述中，我们知道所有的文本单词在生成之前都会先生成它的话题，因此生成的话题序列 $\boldsymbol{z}$ 和文本单词序列 $\boldsymbol{w}$ 是一一对应的。LDA 模型的整体是由观测变量和隐变量组成的联合概率分布，可以表示为:

\[p(\boldsymbol{w,z},\theta, \varphi | \alpha, \beta) = \prod_{k=1}^K p(\varphi_k|\beta) \prod_{m=1}^Mp(\theta|\alpha)\prod_{n=1}^{N_m}p(z_{mn}|\theta_m)p(w_{mn}|z_{mn},\varphi)\]

我们对所有的因变量求积分后得到超参数 $\alpha$ 和 $\beta$ 给定条件下的所有文本生成的概率是:

\[p(\boldsymbol{w} | \alpha, \beta) = \prod_{k=1}^K \int p(\varphi_k|\beta) \left[\prod_{m=1}^M \int p(\theta_m | \alpha) \prod_{n=1}^{N_m} \left[\sum_{l=1}^K p(z_{mn}=l|\theta_m)p(w_{mn}|\varphi_l) \right]d\theta_m \right] d_{\varphi_k}\]

LDA 吉布斯采样算法

LDA模型的学习通常采用收缩的吉布斯抽样，其基本想法是，通过对隐变量 $\theta$ 和 $\varphi$ 积分，得到边缘概率分布 $p(\boldsymbol{w,z}| \alpha, \beta）$ 其中变量 $\boldsymbol{w}$ 是可观测的，变量 $\boldsymbol{z}$ 是不可观测的; 对后验概率分布 $p(\boldsymbol{z} |\boldsymbol{w},\alpha,\beta)$ 进行吉布斯抽样，得到该分布的样本集合；再利用这个样本集合对参数 $\theta$ 和 $\varphi$ 进行估计，最终得到LDA模型 $p(\boldsymbol{w,z},\theta, \varphi | \alpha, \beta)$ 的所有参数估计。

这里省略具体的推导过程直接给出抽样分布的公式:

\[p(\boldsymbol{z} | \boldsymbol{w}, \alpha, \beta) \propto \prod_{k=1}^K \frac{B(n_k + \beta)}{B(\beta)} \cdot \prod_{m=1}^M \frac{B(n_m + \alpha)}{ B(\alpha)}\]

$n_k={ n_{k1}, n_{k2}, \dots, n_{kV}}$，$n_{kv}$ 表示第 $k$ 个话题下，单词 $v$ 出现的次数；其中 $n_m={ n_{m1}, n_{m2}, \dots, n_{mK}}$， $n_{mk}$ 表示第 $m$ 个文本中，话题 $k$ 实际出现的次数;

既然都省略到这种程度了(其实是看证明看的头秃)，那直接给出吉布斯每个分量抽样的条件概率分布表达式:

\[p(z_i| \boldsymbol{z_{-i}}, \boldsymbol{w}, \alpha, \beta) \propto \frac{n_{kv} + \beta_v}{ \sum_{v=1}^V (n_{kv} + \beta_v)} \cdot \frac{n_{mk} + \alpha_k}{\sum_{k=1}^K(n_{mk} + \alpha_k)}\]

注意，这里其中第 $m$ 个文本的第 $n$ 个为止的单词 $w_i$ 是单词集合的第 $v$ 个单词，其话题 $z_i$ 是话题集合的第 $k$ 个话题，$n_{kv}$ 表示第 $k$ 个话题中第 $v$ 个单词的计数减去当前单词的计数，$n_{mk}$ 表示第 $m$ 个话题中第 $k$ 个话题的计数减去当前单词的话题的计数。

当我们抽样结束后，我们就得到了 $\boldsymbol{z}$ 的值，用这些值就可以估计变量 $\theta$ 和 $\varphi$ 啦。对于 $\theta_m$ 我们再来回忆一下它的定义是表示在第 $m$ 个文本中，每个话题 $k$ 出现的概率，并且它的先验服从 $Dir(\theta_m|\alpha)$，我们可以得到它的后验概率，以及其对应的参数估计式:

\[p(\theta_m | \boldsymbol{z_m},\alpha) = \frac{1}{Z_{\theta_m}} \prod_{n=1}^{N_m} p(z_{mn}|\theta_m)p(\theta_m|n_m + \alpha)\] \[\theta_{mk} = \frac{n_{mk} + \alpha_k}{\sum_{k=1}^K (n_{mk} + \alpha_k)}, \quad m=1,2, \dots, M; k=1,2, \dots, K\]

$n_m$ 表示第 $m$ 个文本的不同话题计数。 $\theta_{mk}$ 也就可以看成是它的先验计数与实际它的话题计数之和除以正则项。

类似的，参数 $\varphi$ 也就可以表示为:

\[p(\varphi_k | \boldsymbol{w,z}, \beta) = \frac{1}{Z_{\varphi_k}} \prod_{i=1}^Ip(w_i|\varphi_k)p(\varphi_k|\beta) = Dir(\varphi_k|n_k + \beta)\]

参数估计表达式:

\[\varphi_{kv} = \frac{n_{kv} + \beta_v}{\sum_{v=1}^V(n_{kv}+\beta_v)}, \quad k=1,2,\dots,K; v=1,2, \dots, V\]

其中 $n_k$ 表示第 $k$ 个话题下的单词计数。

吉布斯抽样流程 [2]

输入: 文本单词序列 $\boldsymbol{w = {w_1, \dots, w_m, \dots, w_M}}, \boldsymbol{w_m}=(w_{m1}, \dots, w_{mn}, \dots, w_{mNm})$, 目标概率分布 $p(\boldsymbol{z | w}, \alpha, \beta)$ 超参数 $\alpha$， $\beta$，话题个数 $K$ 输出: 采样的文本话题序列 $\boldsymbol{z={z_1, \dots, z_m, \dots, z_M}}，z_m=(z_{m1}, \dots, z_{mn}, \dots, z_{mN_m})$ 和模型的参数 $\varphi$ 和 $\theta$

初始化所有计数矩阵的元素 $n_{mk}$，$n_{kv}$，计数向量的元素 $n_m$, $n_k$ 的初始值为 0
For $\boldsymbol{w_m} in $ $w_m, m=1,2, \dots, M$
- For $w_{mn}$ in $\boldsymbol{w_m}, n=1,2,\dots, N_m$ (这个文本下的所有单词)
  - 抽样话题 $z_{mn} = z_k \sim Mult(\frac{1}{K})$ (每个话题的概率都为$\frac{1}{K}$)
  - 增加文本-话题计数 $n_{mk} = n_{mk} + 1$
  - 增加文本-话题和计数 $n_m = n_m + 1$
  - 增加话题-单词计数 $n_{kv} = n_{kv} + 1$
  - 增加话题-单词和计数 $n_k = n_k + 1$
循环执行以下操作, 直到进入燃烧期
- For $\boldsymbol{w_m} in $ $w_m, m=1,2, \dots, M$
- For $w_{mn}$ in $\boldsymbol{w_m}, n=1,2,\dots, N_m$ (这个文本下的所有单词)
  - 当前的单词 $w_{mn}$ 是第 $v$ 个单词，话题 $z_{mn}$ 是第 $k$ 个话题；减少计数 $n_{mk} = n_{mk} -1, n_m = n_m-1, n_{kv} = n_{kv} - 1, n_k = n_k-1;$
  - 按照满条件分布进行抽样
  \[p(z_i| \boldsymbol{z_{-i}}, \boldsymbol{w}, \alpha, \beta) \propto \frac{n_{kv} + \beta_v}{ \sum_{v=1}^V (n_{kv} + \beta_v)} \cdot \frac{n_{mk} + \alpha_k}{\sum_{k=1}^K(n_{mk} + \alpha_k)}\]
  (这里可以想象等式的右边就是每个话题的在当前文本$m$和当前单词$v$下面的概率，根据这个概率抽样一个话题) 抽样得到了新的第 $k’$ 个话题，分配给 $z_{mn}$;
  - 增加计数 $n_{mk’} = n_{mk’} + 1, n_m=n_m+1, n_{k’v}=n_{k’v}+1, n_{k’}=n_{k’}+1;$
  - 得到更新的两个计数矩阵 $ N_{K \times V} = [n_{kv}]$ 和 $N_{M \times K} = [n_{mk}]$, 表示后验概率分布 $p(\boldsymbol{z |w},\alpha,\beta)$ 的样本计数；
利用得到的样本计数，计算模型参数:

\[\theta_{mk} = \frac{n_{mk} + \alpha_k}{\sum_{k=1}^K (n_{mk} + \alpha_k)}, \quad m=1,2, \dots, M; k=1,2, \dots, K\] \[\varphi_{kv} = \frac{n_{kv} + \beta_v}{\sum_{v=1}^V(n_{kv}+\beta_v)}, \quad k=1,2,\dots,K; v=1,2, \dots, V\]

简单的代码例子

#-*- coding:utf-8 -*-
import numpy as np
import os

class Lda(object):
    def __init__(self, K, V, alpha, beta):
        self.K = K
        self.V = V
        self.alpha = alpha
        self.beta = beta
        assert len(alpha) == K
        assert len(beta) == V
    
    def learn(self, passages, iterations = 100):
        """
        passages: a list of passage, one passage is a list of words (ids)
        vocabulary: a list of vocabulary (ids)
        """
        # initialize
        M = len(passages)
        V = self.V
    
        self.num_mk = np.zeros((M, self.K))
        self.num_m = np.zeros(M)
        self.num_kv = np.zeros((self.K, V))
        self.num_k = np.zeros(self.K)
        self.topics = []
        for m, passage_m in enumerate(passages):
            cur_topic = []
            for n, word_mn in enumerate(passage_m):
                z_mn = np.random.choice(a=range(self.K), p=[1.0/self.K]*self.K)
                self.num_mk[m][z_mn] += 1
                self.num_m[m] += 1
                self.num_kv[z_mn][word_mn] += 1
                self.num_k[z_mn] += 1
                cur_topic.append(z_mn)
            self.topics.append(cur_topic)
        
        while iterations > 0:
            for m, passage_m in enumerate(passages):
                for n, word_mn in enumerate(passage_m):
                    v = word_mn
                    k = self.topics[m][n]
                    self.num_mk[m][k] -= 1
                    self.num_m[m] -= 1
                    self.num_kv[k][v] -= 1
                    self.num_k[k] -= 1
                    k_plus = self._sample(v, m)
                    self.num_mk[m][k_plus] += 1
                    self.num_m[m] += 1
                    self.num_kv[k_plus][v] += 1
                    self.num_k[k_plus] += 1
                    self.topics[m][n] = k_plus
            iterations -= 1
            if iterations % 5 == 0: 
                print(iterations)
        
        # calculate theta
        temp1 = np.asarray(self.num_mk, dtype=np.float32) + np.asarray(self.alpha)
        self.theta = temp1 / np.sum(temp1, axis=1)[:, np.newaxis]

        # calculate varphi
        temp1 = np.asarray(self.num_kv, dtype=np.float32) + np.asarray(self.beta)
        self.varphi = temp1 / np.sum(temp1, axis=1)[:, np.newaxis]
    
    def _sample(self, v, m):
        p = []
        #normalization_alpha = 0  # alpha可以不用计算，相同的分母项
        for k in range(self.K):
            cur_prob = (self.num_kv[k][v] + self.beta[v]) * (self.num_mk[m][k]+self.alpha[k])
            normalization_beta = np.sum(self.num_kv[k] + self.beta) # 每个话题k下面的word的正则项都不同
            p.append(float(cur_prob)/normalization_beta)  
            #normalization_alpha += self.num_mk[m][k] + self.alpha[k]
        #p = np.array(p) / float(normalization_alpha))

        p = np.array(p, dtype=np.float32)
        normalization = sum(p)
        p = p/normalization
        try:
            k_plus = np.random.choice(a=range(self.K), p = p)
        except:
            print(p)
            
        return k_plus

References

[1] http://ukdataservice.ac.uk/media/307220/presentation4.pdf

[2] 李航. 统计学习方法第二版[J]. 2019.