关于Markov Chain & Monte Carlo

前言：由于近期项目需求，对相关知识做了一下回顾总结。本文仍时以大神[2]的讲解为主线所作的学习笔记，同时参考了很多资料。如有冒犯和疏漏，欢迎批评指正。

Markov Chain & Monte Carlo (MCMC)是推断(Inference)中近似推断中的随机推断。

Monte Carlo Method

Monte Carlo Method是对一类随机方法的特性的概括，即那些“采样越多，越近似最优解”的方法[1]。

Monte Carlo Method是一种基于采样的随机近似方法。推断的任务是求后验概率(posterior probability)，其中为观测变量(observed variable)，为潜在变量(latent variable)。拿到后验后，比较直接的目的是得到其期望，比如函数关于该分布的期望：。但往往关于期望的积分难以求解，所以我们可以通过采样的方式做一种近似，使其从连续的(continuous)变为离散的(discrete)。这就是关于蒙特卡洛方法的思想。

那么问题来了，当后验非常复杂时如何采样？

基于概率分布采样。

要想再后验中采样，一种办法是通过概率密度函数(pdf)求出累积分布函数(cdf)，再加上计算机可以比较容易地在(0,1)之间均匀的生成随机数(cdf纵坐标)，所以可以进一步反推出一些样本x。大体如下丑图。但是问题在于，的过程，对于较复杂的后验来说，难以实现。

拒绝采样(Rejection Sampling)

这里面有一个Proposal Distribution的概念。即假设有这个分布。且有：为常数。为什么要有个呢，因为pdf曲线下面积为1，所以不加限制的话，不会处处大于。那为什么要有的存在呢，因为我们希望有一个简单的分布来代替复杂的，进而通过对进行采样，来判断该采样能否代表。比如，该样本点落入了曲线下则接收。

定义接收率。

方法步骤如下：

  1. 采样：
  2. 利用均匀采样：（为什么要有这个u呢？）
  3. 判断：：接收；：拒绝

重要性采样(Importance Sampling)

不直接对采样，而是对其期望进行采样。和普通的基于后验的采样一样，我们仍然离散化，但这里的离散化不是对的，而是仍然借助了proposal distribution。

注意上式中，。称为Weights。但是这个Weights还是有问题，当二者分布在某点差异大时，有可能会对那些错误的样本分配过高的权重。直观来说：方法的初衷在于，在二者接近时更多的采样。

Markov Chain

时间和状态都是离散的随机过程

齐次(一阶)Markov Chain：当前状态只和上一时刻的状态有关，即

定义状态转移矩阵

我们以概率图的方式考量马氏链。

图中每一个圆圈表示一个对应时刻所观测到的数据，我们不断地观测这些数据地目的在于，希望利用这些数据去找到后验分布，那么如果每一个时刻所观测到的分布都接近于同一个分布或者说都一样，这种情况就是达到了平稳分布。Markov Chain的思路就是希望收敛到平稳分布。

另外有一个平稳分布的充分不必要条件：Detailed Balance。即满足Detailed Balance条件就满足平稳分布条件。注意表示条件概率。

满足Detailed Balance的Markov Chain被称为可翻转的(reversible)。

证明如下：

其中为边缘概率，比较好理解，因为状态转移矩阵的每一行表示该状态转移到其他各状态的概率，所以积分为一。这里是用积分表示的，PRML中是用加和形式表示的，个人认为PRML中更严谨，因为MCMC中对时间和状态的假设都是离散的。这里为了连贯[2]的讲解保持积分形式。

Markov Chain & Monte Carlo

再次明确最终目的：推断后验，进而求其期望。由于直接从中采样困难，所以我们需要间接的去推断。

回顾平稳条件，其中视为，那么转移矩阵如何求？

Metropolis-Hastings(MH算法)

假设我们有一个状态转移矩阵，可称为Proposal Matrix。

有：

因为是proposal的，所以此时左右式并不相等。那么目的就是找到使其相等的状态转移矩阵了。

我们可以通过在左右两端同时乘以一个因子使其两端相等。

简化为：

其中称为接收率。。此时：

所以当接收率这样设计的时候，就达到了Detail Balance条件，从而满足平稳分布。这就是MH算法了。

方法步骤如下：

  1. 采样：
  2. 利用均匀采样：（为什么要有这个u呢？）
  3. 判断：：；：

Gibbs吉布斯采样(特殊的MH采样)

在面对高维分布时，当采样第i维时，固定其它维度(一维一维的采样)。记为。

举例：

  1. t时刻：初始值
  2. t+1时刻：

那么为什么说Gibbs是特殊的MH算法呢？

我们来看中的。对于，对应于Gibbs，显然可以写作，对，可以写成。类似的我们有：

这里面为什么能约掉呢，因为Gibbs在计算第i维度时，其余维度时固定的，也就是和上一个时刻相同，所以就可以约掉了。那么对比下来，Gibbs就是为1的MH算法了。还有就是它没有用proposal的Q了，而是直接哟个了P。

实战

本部分参考：https://kieranrcampbell.github.io/blog/2016/05/15/gibbs-sampling-bayesian-linear-regression.html

这里要估计的模型是：

等价于：

数据：.

目标：求三者的后验，就是方差的倒数。

最大似然：

设定先验：

首先看的更新

有正比关系：

这里面就刚好是上面的似然，另外。

对于正太分布

的话，我们取，且只保留与相关部分，得到：

那么对于来说就是，代入，化简，只保留与相关的项：

展开有：

那么此时的系数就是。的系数为。

综上我们可以定义函数：

def sample_beta_0(y,x,beta_1,tau,mu_0,tau_0):
    N = len(y)
    assert len(x) == N
    precision = tau_0 + tau * N
    mean = tau_0 * mu_0 + tau * np.sum(y - beta_1 * x)
    mean /= precision
    return np.random.normal(mean, 1 / np.sqrt(precision))

然后看更新

和类似，这里省略过程，有：

同样我们定义函数：

def sample_beta_1(y,x,beta_0,tau,mu_1,tau_1):    N = len(y)    assert len(x) == N    precision = tau_1 + tau * np.sum(x * x)    mean = tau_1 * mu_1 + tau * np.sum( (y - beta_0) * x)    mean /= precision    return np.random.normal(mean, 1 / np.sqrt(precision))

再看更新

这里就不是高斯分布了而是分布。

关于分布移步维基百科--伽玛分布。

那么这里概率密度。取对数：。

同上面我们有：

此时概率密度为：

的系数是。的系数为。

那么有：

同样的：

def sample_tau(y, x, beta_0, beta_1, alpha, beta):
    N = len(y)
    alpha_new = alpha + N / 2
    resid = y - beta_0 - beta_1 * x
    beta_new = beta + np.sum(resid * resid) / 2
    return np.random.gamma(alpha_new, 1 / beta_new)

接下来为了进行采样，合成一些数据，以及设定一些超参数设置。

beta_0_true = -1beta_1_true = 2tau_true = 1

N = 50x = np.random.uniform(low = 0, high = 4, size = N)y = np.random.normal(beta_0_true + beta_1_true * x, 1 / np.sqrt(tau_true))

synth_plot = plt.plot(x, y, 'o')plt.xlabel('x')plt.ylabel('y')
# 初始值
init = {'beta_0': 0,
        'beta_1': 0,
        'tau': 2}

# 超参
hypers = {'mu_0': 0,
         'tau_0': 1,
         'mu_1': 0,
         'tau_1': 1,
         'alpha': 2,
         'beta': 1}

def gibbs(y, x, iters, init, hypers):    assert len(y) == len(x)    beta_0 = init['beta_0']    beta_1 = init['beta_1']    tau = init['tau']

    trace = np.zeros((iters, 3)) ## trace to store values of beta_0, beta_1, tau

    for it in range(iters):        beta_0 = sample_beta_0(y, x, beta_1, tau, hypers['mu_0'], hypers['tau_0'])        beta_1 = sample_beta_1(y, x, beta_0, tau, hypers['mu_1'], hypers['tau_1'])        tau = sample_tau(y, x, beta_0, beta_1, hypers['alpha'], hypers['beta'])        trace[it,:] = np.array((beta_0, beta_1, tau))

    trace = pd.DataFrame(trace)    trace.columns = ['beta_0', 'beta_1', 'tau']

    return trace

参考：

https://www.zhihu.com/question/20254139

https://space.bilibili.com/97068901/video

赞 (0)

相关推荐