背景

贝叶斯统计推断是一个非常艺术的东西，他的先验结果（prior knowledge）是一种专家的经验，而贝叶斯公式给出的后验分布是贝叶斯推断的基本工具。每个大学生都应该学过数理统计，其中提及的贝叶斯统计推断都是一些基于简单先验和简单后验的结果。

但当我们一旦碰到复杂的先验（比如高维参数、复杂先验分布），我们用贝叶斯公式得出的后验分布将变得非常复杂，在计算上会非常困难。为此，先人提出了MCMC算法方便我们可以对任何后验分布进行计算或推断。其思想就是其名字：两个MC。

第一个MC: Monte Carlo(蒙特卡洛)。这个简单来说是让我们使用随机数（随机抽样）来解决计算问题。在MCMC中意味着：后验分布作为一个随机样本生成器，我们利用它来生成样本(simulation)，然后通过这些样本对一些感兴趣的计算问题（特征数，预测）进行估计。

第二个MC：Markov Chain(马尔科夫链)。第二个MC是这个方法的关键，因为我们在第一个MC中看到，我们需要利用后验分布生成随机样本，但后验分布太复杂，一些package中根本没有相应的随机数生成函数（如rnorm(),rbinom()）怎么办？答案是我们可以利用Markov Chain的平稳分布这个概念实现对复杂后验分布的抽样。

Markov Chain

为了能顺利阐述MCMC算法，这里就简单地讲一下所涉及到的马尔科夫链概念。因为都是我个人的理解，这里所讲不敢说都是正确的，希望各位童鞋能够独立思考。

1.什么是Markov Chain

随机过程\(\{X_n\},X_n\)的状态空间为有限集或可列集，比如当\(X_n=i\),就称过程在n处于状态i。 定义：如果对任何一列状态\(i_0,i_1,...,i_{n-1},i,j\)，及对任何n≥0，随机过程\(\{X_n\}\)满足Markov性质： \[P\{X_{n+1}=j|X_0=i_0,...X_{n-1}=i_{n-1},X_n=i\}\\=P\{X_{n+1}=j|X_n=i\}\]

转移概率
\(P\{x_{n+1}=j|X_n=i\}\)成为Markov链的一步转移概率\(P_{ij}^{n,n+1}\),当这个概率与n无关,则记之为\(P_{ij}\)

转移概率矩阵P
这个矩阵的元素就是一步转移概率\(P_{ij}\)

结论
一个Markov链可以由它的初始状态以及转移概率矩阵P完全确定。（证明略，自行百度或翻书）

2.n步转移概率

所谓n步转移概率就是从状态i走n步正好到状态j的概率，我们记为\(P_{ij}^{(n)}\)。 利用概率分割的思想，由基础概率论中全概率公式可以得到\[P_{ij}^{(n)}=\Sigma_{k=0}^{\infty}P_{ik}P_{kj}^{(n-1)}\] 写成矩阵形式就是\[P^{(n)}=P\times P^{(n-1)}\] 进一步推广，我们就推出了著名的Chapman-Kolmogorov方程：\[P_{ij}^{(n+m)}=\Sigma_{k=0}^{\infty}P_{ik}^{(n)}P_{kj}^{(m)}\]写成矩阵形式就是\[P^{(m+n)}=P^{(m)}\times P^{(n)}\]

3.Markov链的极限定理和平稳分布

现在我们已经有了n步转移概率的概念，一个很简单的想法就是如果n趋向无穷会怎么样？这个问题也是后面极限分布以及平稳分布的来源，更是MCMC算法中第二个MC的关键。
要回答这个问题首先要掌握几个关键概念，我先列出来，如果不熟悉的可以自行百度或翻书：

  1. 互达性(\(i \leftrightarrow j\))
  2. 周期性(d(i))
  3. 常返与瞬过(\(f_{ii}=1\))
  4. 常返时(\(T_i\))

3.1互达性、周期性

几个重要结论（证明略，自行百度，或者call me）
1. \(i \leftrightarrow j,\Rightarrow d(i)=d(j)\)
2. \(若存在d(i)<\infty,则存在N，对所有n>N,有P_{ii}^{(nd(i))}>0\)
3. \(P_{ji}^{(m)}>0,\Rightarrow 存在N，对所有n>N,有P_{ii}^{(m+nd(i))}>0\)
4. 对于非周期不可约Markov链的转移概率矩阵P，存在N，当n≥N时，\(P^{(n)}>0\)

3.2常返、瞬过、常返时

  1. 引入重要概率\(f_{ij}^{(n)}\)表示从i出发在n步首次到达j的概率，我们约定\(f_{ii}^{(0)}=0\)

  2. \(定义f_{ij}=\Sigma_{n=1}^{\infty}f_{ij}^{(n)}\)表示从i出发最终到j的概率。

  3. \(定义：如果f_{ii}=1,则状态i是常返的，否则是瞬过的\)

  4. \(我们记T_i为首次返回i的步长,且\mu_i=ET_i\)

  5. \(若\mu_i=\infty,则称状态i是零常返，否则正常返\)

3.3 极限定理

是时候回答上面那个问题了！（摘自网络共享PPT）

另外对于一个正常返、非周期的状态(也称为遍历ergodic)，我们有结论:\[对所有i\leftrightarrow j,有limP_{ji}^{(n)}=limP_{ii}^{(n)}=\frac{1}{\mu_i}\]

到这里，我们可以想象了，n步转移概率矩阵最终的极限形式应该是由相同的行向量组成的！！！我们把那个行向量称为极限概率分布。 但是，问题还没完，要计算极限概率，就要根据极限定理，求出很多东西（如\(f_{ii}^{(n)}\)），实际中并不方便。所以，我们要引入平稳分布这个概念。最关键的部分来了！！！！

3.3 平稳分布

什么是平稳分布？它和求极限概率分布有什么关系呢？

定义：Markov链有转移概率矩阵P，如果有一个概率分布\(\{\pi_i \}满足\pi_j =\Sigma_{i=0}^{\infty}\pi_i P_{ij}\),则称为这个Markov链的平稳分布。这个定义用矩阵形式写出来就是π*P=π.

这个定义内容很丰富：如果一个过程的初始状态\(X_0\)有平稳分布π,我们可以知道对所有n，\(X_n\)有相同的分布π。再根据Markov性质可以得到，对任何k，有\(X_n,X_{n+1},...,X_{n+k}\)的联合分布不依赖于n，显然这个过程是严格平稳（这个结论很重要哦）的，平稳分布也由此得名！！

重要定理

若一个Markov链中所有的状态都是遍历的，则对所有i,j有\(limP_{ij}^{(n)}=\pi_j存在，且\pi=\{\pi_j,j≥0\}就是平稳分布\)。反之，拖一个不可约Markov链只存在一个平稳分布，且这个Markov链的所有状态都是遍历的，则这个平稳分布就是该Markov链的极限概率分布。\[limP_{ij}^{(n)}=\pi_j\]

上面这条定理就给出了通过求平稳分布来求Markov链极限概率分布的简单方法。到这里我对第二个MC的解读也就结束了。下一期我们就要开始具体讲讲怎么应用MCMC算法对后验分布进行抽样和推断。

MCMC算法

在前一期的最后3.3一节，我们得到一个重要定理。它讲述了只要有一个Markov Chain是不可约、遍历的，那么它的极限概率分布就可以用它的平稳分布求得。那MCMC算法就是要构建并模拟这么一个Markov Chain使得它是不可约、遍历的而且它的平稳分布就是我们所关心的后验分布。然后，我们就可以以Monte Carlo思想通过这一个Markov Chain的模拟样本代替后验分布做统计分析。

此外，我们还知道，要描述一个Markov Chain只需要给定它的转移概率矩阵和初始值（或者初始分布）就Ok了。所以，MCMC算法的关键就是如何设定这个转移概率矩阵，从中不断抽样，模拟出一个Markov Chain样本。接下来我就给大家介绍两个有名的实现MCMC的算法(参考自维基百科)。

1.Metropolis-Hastings

[Wiki介绍]http://en.wikipedia.org/wiki/Metropolis–Hastings_algorithm

The Metropolis–Hastings algorithm can draw samples from any probability distribution \(g(\theta)\)(这里就是我们所关心的后验分布的密度函数), provided you can compute the value of a function f(x) which is proportional to \(g(\theta)\). The lax requirement that f(x) should be merely proportional to the density, rather than exactly equal to it, makes the Metropolis–Hastings algorithm particularly useful, because calculating the necessary normalization factor is often extremely difficult in practice. (只需要算出后验分布的核函数f(x)就行了哦！)

The Metropolis–Hastings algorithm works by generating a sequence of sample values in such a way that, as more and more sample values are produced, the distribution of values more closely approximates the desired distribution, P(x). These sample values are produced iteratively, with the distribution of the next sample being dependent only on the current sample value (这一串样本就是一个Markov Chain的实现).

Specifically, at each iteration, the algorithm picks a candidate \(\theta^{*}\) for the next sample value based on the current sample value. Then, with some probability P(acceptance probability), the candidate is either accepted (in which case the candidate value is used in the next iteration) or rejected (in which case the candidate value is discarded, and current value is reused in the next iteration).

步骤
1. Choose an arbitrary point \(\theta_0\) to be the first sample.
2. 选择一个概率密度函数(proposal density \(p(\theta^{*}|\theta^{t-1})\))用来生成candidate \(\theta^{*}\)
3. 计算acceptance probability\[P=min\{\frac{g(\theta^{*})p(\theta^{t-1}|\theta^{*})}{g(\theta^{t-1})p(\theta^{*}|\theta^{t-1})},1\}\]
4. 在Markov Chain 的第t步，sample a value \(\theta^t\) and sample \(\theta^{t}=\theta^{*}\) with probability P otherwise \(\theta^t=\theta^{t-1}\)

在MH算法中，对proposal density的不同取法使得MH自身在各种软件包中有不同的命令名。比如R中就有rwmetrop,indepmetrop

2.Gibbs Sampling

另外还有一个更加有名的算法叫吉布斯抽样，它其实是Metropolis–Hastings算法的一个特例,下面我们考虑多维参数\(\mathbf{X} = (x_1, \dots, x_n)\)。

思想

The point of Gibbs sampling is that given a multivariate distribution it is simpler to sample from a conditional distribution than to marginalize by integrating over a joint distribution. Suppose we want to obtain k samples of \(\mathbf{X} = (x_1, \dots, x_n)\) from a joint distribution \(\left.p(x_1, \dots, x_n)\right.\). Denote the ith sample by \(\mathbf{X}^{(i)} = (x_1^{(i)}, \dots, x_n^{(i)})\). We proceed as follows:

We begin with some initial value \(\mathbf{X}^{(0)}\). To get the next sample (call it the i+1th sample for generality) we sample each component variable \(x_j^{(i+1)}\) from the distribution of that variable conditioned on all other variables, making use of the most recent values and updating the variable with its new value as soon as it has been sampled. This requires updating each of the component variables in turn. 比如，If we are up to the jth component\(x_j\), we update it according to the distribution specified by \(p(x_j|x_1^{(i+1)},\dots,x_{j-1}^{(i+1)},x_{j+1}^{(i)},\dots,x_n^{(i)})\). Note that we use the value that the j+1th component had in the ith sample not the i+1th.此外，在实际操作中我们可以对抽样的结果再添加一个噪音Z上去\(x_j^{*}=x_j+c_j\times Z\)

Repeat the above step k times.

如果在利用\(p(x_j|x_1^{(i+1)},\dots,x_{j-1}^{(i+1)},x_{j+1}^{(i)},\dots,x_n^{(i)})\)对参数向量中的未知参数\(x_j\)抽样时感到困难，可以利用之前提到的Metropolis-Hastings算法代替。这种办法，我们称之为“a Metropolis within Gibbs”

有以下结论:

  1. The samples approximate the joint distribution of all variables.
  2. The marginal distribution of any subset of variables can be approximated by simply considering the samples for that subset of variables, ignoring the rest.
  3. The expected value of any variable can be approximated by averaging over all the samples.

注意:

  1. The initial values of the variables can be determined randomly or by some other algorithm such as expectation-maximization.
  2. It is not actually necessary to determine an initial value for the first variable sampled.
  3. It is common to ignore some number of samples at the beginning (the so-called burn-in period), and then consider only every nth sample when averaging values to compute an expectation. 这么做的原因是 (1) successive samples are not independent of each other but form a Markov chain with some amount of correlation; (2) the stationary distribution of the Markov chain is the desired joint distribution over the variables, but it may take a while for that stationary distribution to be reached. Sometimes, algorithms can be used to determine the amount of autocorrelation between samples and the value of n (the period between samples that are actually used) computed from this, but in practice there is a fair amount of “black magic” involved.

3.MCMC推断

以上，我们已经获得了平稳分布是所关心的后验分布的Markov Chain的一个实现样本。接下来，我们就可以通过trace plot剔除一些因为进化需要生成的开始阶段的样本(burn-in period sample)。然后，我们就可以像对待平稳过程样本一样对待这个序列样本做推断，我就不多说喽！！！

讲完啦，今晚欧冠决赛我看好巴萨，你呢？