由LDA模型，追溯贝叶斯学派的概率情结和采样方法的来龙去脉

LDA

LDA简述

隐狄利克雷分布（Latent Dirichlet Allocation,缩写LDA）是2003年由David M.Blei,Andrew Y.Ng ,Michael I.Jordan 提出的，是一种词袋模型（也是lDA的天然缺陷，不像lstm,transform一类模型，考虑前面顺序，但在深度学习未兴起之前，LDA还是以其的复杂性和不错的表征能力吸引了很多关注），它认为文档是一组词构成的集合，词与词之间无序，一篇文档可以包含多个主题，文档中的每个词都是由某个主题生成的，LDA的结果是能给出每个文档下每个主题的概率分布，同时给出每个主题上词的概率分布，LDA是一种无监督学习，能在文档主题识别，文本分类，文本相似度计算和文章相似度推荐方面有所应用。
这便是LDA的idea,下面会从贝叶斯公式，讲到plsa,再讲到LDA，后面还会讲到LDA求解用到的吉布斯采样和EM算法（解决含隐变量问题的常规标准操作算法），讲好吉布斯采样就又会需要从马尔可夫链讲起，毕竟吉布斯采样属于马尔可夫蒙特卡洛采样方法（mcmc）的范畴内，会讲明白吉布斯采样究竟在干一件什么事，目的是什么，已经吉布斯采样为什么可以这么做，由什么原理保证的（吉布斯采样的推导）

贝叶斯

先从贝叶斯公式出发：（为什么是贝叶斯，因为在LDA的发展路径上，都是贝叶斯学派在搞事情，纵观LDA的诞生，就是一个不断把一个变量变成估计一个随机变量的过程，所有的过程被看成是概率随机的 ）

$p(\theta |X)=\frac{p(\theta ,X)}{p(X)}=\frac{p(X|\theta )\ast p(\theta )}{p(X)}$
即就是我们长说的，贝叶斯公式的物理意义：
$posterior=\frac{likelihood\ast prior}{evidence}$

$后验概率=\frac{似然概率\ast 先验概率}{规范化的证据}$
这就是贝叶斯公式的密度函数的形式，即在给定样本X的条件下，$\theta$的条件分布$p(\theta |X)$就被称为$\theta$的后验分布。它是集中了总体（公式的分母），样本（公式的分子第一项），和先验（公式的分子第二项）三种信息中有关$\theta$的一切信息，而又是排除一切与$\theta$无关的信息之后得到的结果，故基于后验分布$p(\theta |X)$对$\theta$进行统计推断是更合理的，用人话来说就是：

1. 如果$\theta$代表我们的分类变量，例如取值0和1，而$X$代表样本，那么先验概率$p(\theta )$是所有样本整体分类的概率，即我们闭着眼睛随机拿一个样本，不观察它，它是0或1的概率，很简单，抛硬币，押大小，都会根据这个先验概率来猜，一般我们都会猜一个先验概率大的结果
2. 后验概率，$p(\theta |X)$是我们会观察当前样本，再去分析它究竟属于哪一个分类，这不就是我们一般工作中要解决的问题，分类
3. 所有所谓基于分布和概率的方法都是在贝叶斯框架下估计后验概率的分布，从而得到我们判断的分数，或者叫概率。

贝叶斯公式的经典之处在于将后验概率的问题，转化为了先验概率（或者说先验概率*似然概率）的问题，显然，一般情况下，先验概率和$X$的分布，我们都可以根据样本估计得到，问题是，$likelihood$怎么求，？在各类先验已知的情况下，似然的大小就决定了分类结果，因为分母$p(X)$是常数，对于比较概率没有意义，基于这样的方法，我们可以统一的将他们归纳为贝叶斯框架下的分类方法，简称为“贝叶斯分类”，这是一种从贝叶斯公式出发，基于统计学习的方法，为什么又是基于统计学习呢？因为最大似然估计里的似然计算方法其实是利用统计函数进行模糊匹配的，我们假！设！似然的$p(X|\theta )$分布是符合正态分布的，这个正态分布的均值估计和方差估计都是从现有的样本进行统计得出的，而最简单最经典的估计方法就是最大似然概率方法。当不考虑先验概率或或者默认先验概率相等下，直接来对
$\theta =0$和$\theta =1$的似然分数进行比较，而比较的这个东西，就是所谓的似然比，likelihood rate（LR），$LR=\frac{p(X|{\theta }_1)}{p(X|{\theta }_2)}>or<1$
如果在ml的基础上，进一步将先验考虑进来，这样的方法就被称为MAP，最大后验概率方法，注意，引入先验概率后的方法被称为最大后验概率方法（单从名字上看，有点怪，老外对这个算法的起名真的不太友好，一点都不见名知义，但本质上没毛病的）。。。因为$p(X)$对比较概率没有意义，这时候的比较，就相当于计算不同类的后验概率$p(\theta |X)$：
$LR=\frac{p({\theta }_1)p(X|{\theta }_1)}{p({\theta }_2)p(X|{\theta }_2)}>or<1$
如果在ml的基础上，引入误判代价，就是所谓的最小风险贝叶斯方法，此处不再详细说了，贝叶斯框架是分类问题最经典的解决方法，这是一个特别严格而且不讲道理的假设，但是它又能够覆盖大多数的情况。

LDA的前菜plsa

回到LDA的问题：从概率分布的角度来解决文章主题的分析问题，这类问题解决最简单的方式是一元模型的方法，unigram model，其主要思想是，上帝按如下的规则产生一个文本：

1. 上帝有一个骰子，这个骰子有V个面，每个面上有一个单词，各个面的概率不同，（v代表vocabulary大小）
2. 上帝开始抛骰子，每抛出一个面就产生一个单词，如果一个文档有n个词，上帝需要独立抛n次

骰子的各个面的概率是：$\vec{p}=(p_1,p_2,p_3,,,p_V)$,$p_i$代表产生每个词的概率，那么对于一篇文档$\vec{w}=(w_1,w_2,w_3,,,w_n)$,$w_i$代表 文档中第$i$个词，生成文档的概率是：
$p(\vec{w})=p(w_1,w_2,w_3,,,w_n)={\prod }_{j=1}^{n}p(w_j)$
我们的目的就是获取上帝玩游戏的信息，即$\vec{p}$的分布，可以通过最大似然概率去求解这个值。

进一步把这个一元模型延伸到贝叶斯一元模型unigram model：
在贝叶斯派看来，一切参数都是随机变量，所以模型中的$\vec{p}$不能是唯一固定的，而是需要服从一个分布，也就是说，上帝不应该只有一个固定的骰子，还是应该有一个盒子，这个盒子中装满了各种骰子，上帝先从盒子里挑出一个骰子，然后再用这个骰子来玩游戏，上帝的这个盒子里，有多种骰子，每种骰子的数量不一样，有些骰子数量多，有些骰子数量少，所以从概率分布的角度看，盒子里的骰子$\vec{p}$，服从一个概率分布$p(\vec{p})$,这个分布是参数$\vec{p}$的先验分布，先验分布$p(\vec{p})$,可以有很多种选择，最好的选择是狄利克雷分布（Dirichlet ），终于出现了狄利克雷分布，（至于为什么最好，有解释说狄利克雷分布是多项式分布的共轭分布，好吧又是一个坑），然后对$p(\vec{p})$应用贝叶斯估计，即似然概率*先验概率的套路，可以得到语料的产生概率。

发展到这里，就有了一个问题，一元模型unigram model没有考虑文档有多个主题的情况，或者说一元模型的本质是把文档看成是只有一个主题的情况下产生的，一般情况下，一篇文档可以看成是由多个主题（Topic）组成的，文档中的每个词又被认为是从由一个固定的主题（Topic）生成的，这就衍生出了LDA发展路径上的另一个重要算法，plsa,plsad的游戏规则是：
1.上帝有两中类型的骰子，一类是doc-topic,这类骰子有k个面，每个面都是一个topic，（这里假设了doc-topic这类骰子只有一个），另一类是topic-word,这个是上帝之前经常玩的骰子
2.上帝的topic-word有k个，对应了doc-topic骰子有k个面。
3.生成每篇文档前，上帝都先为这个文档，制造一个“特定的”骰子，然后重复下面的动作，
3.1 抛这个doc-topic骰子，得到编号为z的topic,
3.2 从k个topic-word骰子中，选出第z个骰子，并抛这个骰子，得到一个面，即产生一个词

PLSA模型中doc-topic 骰子和 topic-word骰子 的每个面的概率都是固定值，所以所以点估计，但是plsa模型中即包含观测变量$d_i,w_j$,能观测到的每篇文档，以及文章中的每个词,又多出了隐变量$z_k$,每个文档是观测不到其中包含的主题，因为隐变量$z_k$的存在，就不能简单的采用极大似然函数的方法求解模型参数了，但对隐变量$z_k$的模型参数求解，也有固定套路，那就是EM算法，这个坑以后再填

于是，我们用EM算法，来估计$p(\vec{p})$中对模型参数，
$D=(d_1,d_2,...,dN)$表示在语料中有N篇文档
$W=(w_1,w_2,...,wM)$表示在语料中有M种词
$n_{ij}(d_i,w_j)$表示在文档$d_i$中$w_j$出现对次数
$Z=(z_1,z_2,...,z_K)$表示$K$个主题，根据3.1的玩法，每篇文档会有多主题

$p(d_i)$表示文档$d_i$的概率
$p(z_k|d_i)$给定文档$d_i$下，主题$z_k$出现的概率，对应3.1
$p(w_j|z_k)$给定主题$z_k$下，单词$w_j$出现的概率，对应3.2

这里一个核心的idea，就是文档和词都是可以直接观测的，而主题是一个隐变量，无法直接观测到，即我们能看到的是“doc-word”（现象）,而本质是“doc-topic-word",即：
$p(d_i,w_j)={\sum }_{k=1}^{K}p(d_i)p(z_k|d_i)p(w_j|z_k)$

具体解法需要结合em算法来讲，大致过程是，固定$p(z_k|d_i)$的分布，求最优的$p(w_j|z_k)$，然后固定$p(w_j|z_k)$，重新计算$p(z_k|d_i)$的分布，并再一次固定新的$p(z_k|d_i)$的分布，再次求解最优的$p(w_j|z_k)$，按此套路迭代，直至收敛，注意，EM算法的收敛性可严格证明。
讲到这里，其实并没有提到吉布斯采样，哈哈，其实plsa还算比较简单的机器学习问题，无非多出来一个隐变量主题的概念，然后有固定套路em算法来解决它，不算很难理解。

下面将进入真正的LDA算法

LDA算法

发展到这里，贝叶斯学派就不同意了，如同从一元模型到plsa的改进，贝叶斯学派对plsa也想作出类似的改进，贝叶斯学派开始认为：为什么上帝只有一个doc-topic的骰子呢，为什么上帝只有固定k个topic-word个骰子呢？于是LDA算法正式出现，
在plsa算法的基础上，玩法进一步升级，复杂
1. 上帝有两个盒子的骰子，一个盒子里装有用于doc-topic的多个骰子，这类骰子每个都有 K个面，每个面对应一个topic,另一个盒子里装有用于topic-word的骰子
2. 上帝随机地从第二个骰子中独立的选择K个骰子，编号为1到K，（注意，上帝不再是一来就有固定的K个骰子，而是先从第二类盒子中挑出了K个骰子，这是与plsa的区别1）
3. 在生成每篇语料之前，上帝都先从第一个盒子中随机抽取了一个doc-topic的骰子（注意，不再是定制的一个了，这是与plsa的区别2），然后重复抛这个骰子，为文档中每一个词都生成了一个topic，（注意，这句话重点理解）
4. 对于语料中每篇文档的每一个编号为z的主题，从K个topic-word骰子中选出第z个骰子，并抛这个骰子，得到一个面，即产生一个词

说到这里，不知道有没有乱，其实，总结一下，LDA和plsa的区别就是：
1. doc-topic的骰子，从一个变成了从盒子中随机拿一个
2. topic-word的骰子，从一组固定的K个变成了从盒子中随机独立地拿K个

由于topic-word从原来的固定一组，变成来随机变量，
${\vartheta }_i=p(z|d_i)$表示在给定文档$d_i$下，主题z出现的概率分布，注意，由从plsa中的概率变成了概率分布
${\phi }_k=p(w|z_k)$表示在给定主题$z_k$下，单词w出现的概率分布，注意，也由从plsa中的概率变成了概率分布

LDA中我们会假设$\vartheta =p(z|d_i)$，${\phi }_k$都是一个狄利克雷分布（Dirichlet ），第二次提狄利克雷分布。
$\vec{\alpha }$：是$\vartheta$先验分布的超参，需要经过这个Dirichlet分布采样得到$\vartheta$，采样方法：吉布斯采样
$\vec{\beta }$：是$\phi$先验分布的超参，也需要经过这个Dirichlet分布采样得到$\phi$，采样方法也是吉布斯采样

所以，在先验条件下，产生主题产生词的联合概率分布为：
$p(\vec{w},\vec{z},|\vec{\alpha },\vec{\beta })=p(\vec{w}|\vec{z},\vec{\beta }).p(\vec{z}|\vec{\alpha })$

根据这个公式，可以推导出（很复杂，后续填坑）参数的贝叶斯公式，再结合吉布斯采样，得到狄利克雷分布样本，从而得到最后的贝叶斯参数估计。

以上就是LDA的整体思想idea和诞生过程，最后我们回到算法的名字，隐狄利克雷分布（Latent Dirichlet Allocation,缩写LDA），再说两点：
1，上帝从两个盒子里取骰子是基于Dirichlet分布的
2，在词生成的马尔可夫过程中（吉布斯采样是基于马尔可夫过程的），主题topic是一个隐变量Latent

填坑1，先填马尔可夫过程和马尔可夫蒙特卡洛采样（MCMC）的坑

马尔可夫链和MCMC

马尔可夫链
马尔可夫链的意义在于，假设有一个连续的状态或过程，下一个状态或过程的发生概率只依赖于前一个状态，与前前一个，以及更久之前的状态以没有关系了，从公式上可以表述成，
$P(X_{t+1}=x|X_t,X_{t-1},...,X_0)=P(X_{t+1}=x|X_t)$

1.假设可能存在3种状态相互转换，我们定义矩阵P是一个状态转移矩阵，其中P(i,j)的值是p(j|i)表示从状态i到状态j的概率
2.假设存在一个初始状态的值${\pi }_0=({\pi }_{01},{\pi }_{02},{\pi }_{03})$,其中，${\pi }_{0i}$表示 在初始时刻，3个状态当前的概率，那么${\pi }_1={\pi }_{0}P$,${\pi }_2={\pi }_{1}P={\pi }_0{P}^2$,…${\pi }_n={\pi }_0{P}^n$
3.可以改变初始状态的值，可以发现一个现象，经过足够多的迭代，当n足够大的时候，总是收敛到同一个概率分布$\pi$，而这个收敛的现象与初始状态的概率分布无关，而是由转移矩阵P决定的。
4.当n足够大以后，$P^n$的每一行都稳定的收敛到了$\pi$。

总结马尔可夫的收敛现象，有马氏链定理：
如果一个非周期的马尔可夫链具有转移矩阵$P$,且它的任意两个状态都是连通的，（从任意一个状态可以通过有限步到达其他任意状态，不会出现条件概率一致为0导致不可达的情况，或者跳不出去的黑洞现象），那么$li{m}_{n->\infty }{P}_{ij}^n$存在且与i无关，并且：

1. $li{m}_{n->\infty }{P}_{ij}^n={\pi }_j$,即$P^n$的每一行都稳定的收敛到同一个$\pi$,$li{m}_{n->\infty }{P}_{ij}^n=\left[\begin{array}{ccccc}\pi (1)& \pi (2)& ...& \pi (j)& ...\\ \pi (1)& \pi (2)& ...& \pi (j)& ...\\ ...& ...& ...& ...& ...\\ \pi (1)& \pi (2)& ...& \pi (j)& ...\end{array}\right]$

2.$\pi (j)={\sum }_{i=0}^{\infty }\pi (i)P_{ij}$ 这个公式解释其意义，当前时刻某一个状态j的概率等于所有状态转移到此状态的概率和
3.$\pi$是方差.$\pi P=\pi$ 的唯一非负解，且${\sum }_{j=1}^{\infty }\pi (j)=1$ 称为马氏链的平稳分布，这也是构造P的方法吗？

如果从一个具体的初始样本$x_0$开始，沿着马氏链按照概率转移矩阵做跳转，那么可以得到一个转移序列，$x_0,x_1,...,x_n,x_{n+1},...$,由于马氏链的收敛行为，$x_n,x_{n+1}，...$都将是平稳分布$\pi$的样本，注意，这里$x_0$是样本，不是一个初始状态
这个马氏链的收敛定理是非常重要的，所有的马尔可夫-蒙特卡洛方法（MCMC）都是以这个定理作为理论基础的，并利用马尔可夫链的平稳分布这个概念。

MCMC

首先我们要明确MCMC是干什么的

一句话总结：MCMC是用来生成样本的。在实际工程中，我们需要利用后验分布生成随机样本，但是后验分布太复杂，简单的均匀分布，高斯正态分布有随机生成函数，基础的package中没有一些服从复杂分布的随机生成函数(只有rnorm(),rbinom())怎么办？答案是我们可以利用马尔可夫链的平稳分布这个概念来实现对复杂后验分布的抽样（样本生成）。这就是MCMC方法。
可以看到LDA中很多步骤 是需要按狄利克雷分布进行“样本”采样的（在这之后才能做参数估计），这就是MCMC方法在LDA中的地位，用于解决复杂分布的采样问题，或者说复杂分布样本生成问题。

MCMC全名：马尔可夫链-蒙特卡洛方法，把马尔可夫（Markov）过程引入到Monte Carlo模拟中，实现抽样分布随着模拟的进行而改变的动态模拟。MCMC的提出是基于下面的问题：对于概率分布p(x)，我们希望得到由p(x)生成的一系列样本，借鉴马氏链定理，如果我们可以构造一个马氏链，它符合：
1. 转移矩阵为$P$的非周期马氏链
2. 这个马氏链是收敛的，且它的平稳分布就是我们的目标分布$p(x)$

那么很显然，只要我们从任何一个初始状态来出发，沿着马氏链不停的转移，如果这个马氏链在第n步收敛，那么我们就可以得到服从$p(x)$的样本$x_n,x_{n+1},...$

所以下一个问题是如何构造转移矩阵P,这需要基于下面的细致平稳条件：
如果转移矩阵P和分布 $\pi (x)$ 满足：
$\pi (i)P_{ij}=\pi (j)P_{ji}$
则，$\pi (x)$将会是马氏链的平稳分布，同时P也将是这个分布的转移矩阵。
证明：
${\sum }_{i=1}^{\infty }{\pi }_i{P}_{ij}={\sum }_{j=1}^{\infty }{\pi }_j{P}_{ji}={\pi }_j{\sum }_{i=1}^{\infty }{P}_i={\pi }_j->\pi P=\pi$

根据$\pi (i)P_{ij}=\pi (j)P_{ji}$这个细致平稳条件，我们在MCMC中的思路就是寻找转移矩阵P,只要P使得概率分布$p(x)$满足细致平稳条件，我们就可以想要的东西：服从$p(x)$分布的一系列样本

假设有一个转移矩阵为Q的马氏链，（这个Q并不是我们最终要的转移矩阵），显然一般都有，　$p(i)Q(i,j)\ne p(j)Q(j,i)$
即细致平稳条件不成立，$p(x)$不是转移矩阵Q的马氏链的平稳分布。不要紧，聪明的我们可以加上一个小技巧，使得平稳条件成立，就让引入一个$\alpha (i,j)$,它真正的名字一般称为接受率，这样的接受率使得：
$p(i)Q(i,j)\alpha (i,j)=p(j)Q(j,i)\alpha (j,i)$

问题又来了，$\alpha (i,j)$怎么取，聪明的我们又看到了对称性，：
$\alpha (i,j)=p(j)Q(j,i)$
$\alpha (j,i)=p(i)Q(i,j)$

接着我们定义表达式：
$Q^{\prime }(i,j)=Q(i,j)\alpha (i,j)$
$Q^{\prime }(j,i)=Q(j,i)\alpha (j,i)$
这样，则：
$p(i)Q^{\prime }(i,j)=p(j)Q^{\prime }(j,i)$

由此，我们好像原来一个具有转移矩阵是Q的马氏链，改造成了转移矩阵是$Q^{\prime }$，且满足细致平稳条件，并且平稳分布是$p(x)$的收敛马氏链。

我们可以得到一个非常好的结论，转移矩阵$Q^{\prime }$可以通过一个马氏转移矩阵Q乘以$\alpha (i,j)$得到，$\alpha (i,j)$一般称为接受率，其取值范围是[0,1],可以理解为一个概率值，看一下公式：
$p(i)Q^{\prime }(i,j)=p(i)Q(i,j)\alpha (i,j)$

原来转移矩阵为Q的马氏链上，从状态i,以$Q(i,j)$的概率转移到状态j的时候，我们以但一定的概率$alpha(i,j)$接受了这个转移（注，这个就是接受-拒绝采样方法里的目标分布和已知分布的比$\frac{\pi (x)}{Mq(x)}$,M是常数，说明了目标分布$\pi (x)$和已知分布$q(x)$的相似性）。

和接受-拒绝采样方法相同的思想，以一个常见的分布通过一定的接受-拒绝概率来得到一个不常见的分布，这里以一个常见的 马氏链的转移矩阵Q，通过一定的接受-拒绝概率得到新的马氏链状态转移矩阵$Q^{\prime }$，（备注，接受-拒绝采样方法可以证明，通过接受-拒绝方法，可以通过一个简单的已知分布得到服从目标分布的样本，后续会增加这部分的说明）。

MCMC采样方法的伪代码 表述如下：
1.初始化马氏链的初始状态$X_0=x_0$
2.for i=0,1,2…,循环一下过程进行采样：
2.1 第t个时刻的马氏状态$X_t=x_t$,基于转移矩阵，采样下一个状态 $y\sim Q(y|x_t)$,
2.2从均匀分布采样 $u\sim Uniform[0,1]$
2,3 如果$u<=\alpha (x_t,y)=p(y)Q(y,x_t)$,则接受转移$x_t->y$,即$X_{t+1}=y$
2.4否则不接受转移，$X_{t+1}=X_t$

有一点请注意，上述流程中的几个点，$y$是我们通过当前转移矩阵$Q$采样得到的样本，$p(y)$是我们通过目标平稳分布得到的样本在$p$上的概率.

个人简单理解，显然$\alpha (x_t,y)=p(y)Q(y,x_t)$越大，转移预热更容易被接受，也就是说，当$p(y)$较大的时候，得到的是平稳分布容易采样大概率得到的样本，这部分样本最容易代表$p(x)$的分布，当$Q(y,x_t)$较大时，表示原马氏链的跳转概率很大，这部分跳转产生的样本有助于式马氏链收敛。综合两部分来看，采样得到的样本是要在符合平稳分布和现有状态之间的均衡，初始状态随机，在状态的不停跳转过程中，会慢慢迁移到这样一个状态，当前状态是一种跳转到$p(y)$采样有较大概率能得到的状态的状态。比较拗口。

现在介绍的MCMC 还是有风险的，万一初始状态不合理，接受率一直很低，跳转概率很低，或者干脆跳不出来，进行进入了黑洞一样，怎么办，那就该介绍M-H采样的故事了，以及吉布斯采样的故事了。
未完待续。。。