最近我们被客户要求撰写关于吉布斯采样的研究报告。

由Kaizong Ye，Liao Bao撰写

我将得出block的Gibbs采样器所需的条件后验分布。

可下载资源

完整程序、数据和文档（word）

贝叶斯模型

假设我们有一个样本量的主题。贝叶斯多元回归假设该向量是从多元正态分布中提取的，通过使用恒等矩阵，我们假设独立的观察结果。

常用于DBM和DBN，吉布斯采样主要用在像LDA和其它模型参数的推断上。

要完成Gibbs抽样，需要知道条件概率。也就是说，gibbs采样是通过条件分布采样模拟联合分布，再通过模拟的联合分布直接推导出条件分布，以此循环。

概念解释

吉布斯采样是特殊的Metropolis-Hastings算法，会用到马尔科夫链。
具体地说，
MCMC：Markov链通过转移概率矩阵可以收敛到稳定的概率分布。这意味着MCMC可以借助Markov链的平稳分布特性模拟高维概率分布；当Markov链经过burn-in阶段，消除初始参数的影响，到达平稳状态后，每一次状态转移都可以生成待模拟分布的一个样本。

Gibbs抽样是MCMC的一个特例，它交替的固定某一维度，然后通过其他维度的值来抽样该维度的值，注意，gibbs采样只对z是高维（2维以上）（Gibbs sampling is applicable in situations where Z has at least two dimensions）情况有效。

吉布斯采样的通俗解释

Gibbs Sampling就是以一定的概率分布，看发生什么事件。

例子

甲只能E：吃饭、学习、打球，
时间；T：上午、下午、晚上，
天气；W：晴朗、刮风、下雨。
现在要一个sample，这个sample可以是：打球+下午+晴朗。

问题是我们不知道p(E,T,W)，或者说，不知道三件事的联合分布joint distribution。当然，如果知道的话，就没有必要用gibbs sampling了。但是，我们知道三件事的conditional distribution。也就是说，p(E|T,W),p(T|E,W),p(W|E,T)。现在要做的就是通过这三个已知的条件分布，再用gibbs sampling的方法，得到联合分布。

具体方法

首先随便初始化一个组合,i.e. 学习+晚上+刮风，
然后依条件概率改变其中的一个变量。
具体说，假设我们知道晚上+刮风，我们给E生成一个变量，比如，学习-》吃饭。我们再依条件概率改下一个变量，根据学习+刮风，把晚上变成上午。类似地，把刮风变成刮风（当然可以变成相同的变量）。这样学习+晚上+刮风-》吃饭+上午+刮风。
同样的方法，得到一个序列，每个单元包含三个变量，也就是一个马尔可夫链。然后跳过初始的一定数量的单元（比如100个），然后隔一定的数量取一个单元（比如隔20个取1个）。这样sample到的单元，是逼近联合分布的。

二维吉布斯采样算法

吉布斯采样算法中右边的条件概率我们是知道的，例如你要采样的是二维高斯分布，那么固定xt后就是二维高斯分布固定xt后的一维高斯分布，且每次采样的坐标不同，这样这个一维高斯分布概率密度函数也就不一样了。

到目前为止，这与环境中看到的多元正态回归相同。则将概率最大化可得出以下解：

贝叶斯模型是通过指定为一个先验分布得到。在此示例中，我将在以下情况下使用先验值

视频

贝叶斯推断线性回归与R语言预测工人工资数据

探索见解 ➜

去bilibili观看 ➜

探索更多视频 ➜

block Gibbs

在对采样器进行编码之前，我们需要导出Gibbs采样器的每个参数的后验条件分布。

条件后验取更多的线性代数。

这是一个非常漂亮和直观的结果。条件后验的协方差矩阵是协方差矩阵的频繁估计，

还要注意，条件后验是一个多元分布，因为它是一个向量。因此，在Gibbs采样器的每次迭代中，我们从后验画出一个完整的矢量。

模拟

我模拟的结果向量。

运行 Gibbs采样器会生成对真实系数和方差参数的估计。运行了500,000次迭代。修整周期为100,000次，修整了10次迭代。

最受欢迎的见解

1.matlab使用贝叶斯优化的深度学习

2.matlab贝叶斯隐马尔可夫hmm模型实现

3.R语言Gibbs抽样的贝叶斯简单线性回归仿真

4.R语言中的block Gibbs吉布斯采样贝叶斯多元线性回归

5.R语言中的Stan概率编程MCMC采样的贝叶斯模型

6.Python用PyMC3实现贝叶斯线性回归模型

7.R语言使用贝叶斯层次模型进行空间数据分析

8.R语言随机搜索变量选择SSVS估计贝叶斯向量自回归（BVAR）模型

9.matlab贝叶斯隐马尔可夫hmm模型实现

以下是MCMC链的图，其中真实值用红线表示。


# 计算后验摘要统计信息（未在其余代码中使用的统计信息）
post_sum_stats<-post_dist %>%
  group_by(param) %>%
  summarise(median=median(draw),
            lwr=quantile(draw,.025),
            upr=quantile(draw,.975)) %>%
  mutate(true_vals=c(tb,tphi))

# 合并汇总统计信息
post_dist <- post_dist %>%
  left_join(post_sum_stats, by='param')

# 绘制MCMC链
ggplot(post_dist,aes(x=iter,y=draw)) +
  geom_line() +
  geom_hline(aes(yintercept=true_vals, col='red'), show.legend=FALSE)+
  facet_grid(param ~ .,scale='free_y',switch = 'y') +
  theme_bw() + 
  xlab('Gibbs Sample Iteration') + ylab('MCMC Chains') + 
  ggtitle('Gibbs Sampler MCMC Chains by Parameter')


ggplot(post_dist,aes(x=draw)) +
  geom_histogram(aes(x=draw),bins=50) +
  geom_vline(aes(xintercept = true_vals,col='red'), show.legend = FALSE) +
  facet_grid(. ~ param, scale='free_x',switch = 'y') +
  theme_bw() + 
  xlab('Posterior Distributions') + ylab('Count') + 
  ggtitle('Posterior Distributions of Parameters (true values in red)')

这是修整后参数的后验分布：