如果您可以写出模型的似然函数，则 Metropolis-Hastings算法可以负责其余部分（即MCMC ）。

由Kaizong Ye，Liao Bao撰写

我写了r代码来简化对任意模型的后验分布的估计。具体如下：

我写了r代码来简化对任意模型的后验分布的估计。具体如下：

可下载资源

 
 
a<-pars[1]      #截距
 
b<-pars[2]      #斜率
 
sd_e<-pars[3]   #残差
 
if(sd_e<=0){return(NaN)}
 
 
 
log_likelihood<-sum( dnorm(data[,2],pred,sd_e, log=TRUE) )

1.1 基本概念

从一个概率分布（目标分布 $P(x)$ ）中得到随机样本序列
这个序列可以用于：a) 近似估计分布 $P(x)$ ; b) 计算积分(期望)
用于高维分布取样
缺陷: MCMC固有缺点, 样本自相关性

1.2 优势

可以从任意的概率分布 $P(x)$ 中取样，只要满足条件：函数 $f(x)$ 成比例于 $P(x)$ 的密度。
更宽松的要求： $f(x)$ 仅需要与 $P(x)$ 的密度成比例。

1.3 要点

生成样本值序列；样本值产生得越多，这些值的分布就越近似于 $P(x)$
迭代产生样本值：下一个样本的分布仅仅取决于当前样本值（马尔科夫链特性）
接受/拒绝概率：接受计算出的值为下一个样本值/拒绝并重复使用当前样本值；基于 $P(x)$ ，接受概率通过比较 $f(x_t)$ （当前值）和 $f(x')$ （备选值）得出

1.4 Metropolis Algorithm (对称分布)

Input: $f(x)$ ，与目标分布 $P(x)$ 成比例的函数

1. 初始化：

$x_0$ : 取任意点作为第一个样本值
$g(x|y)$ : 任意的概率密度函数；给定上一个样本值 $y$ ，计算出下一个样本值 $x$ ；一般我们设定 $g(x|y)$ 是对称的，即 $g(x|y)=g(y|x)$ ；我们也会令 $g(x|y)$ 服从以 $y$ 为中心的 $Gaussian$ 分布，即越靠近 $y$ ，概率越大
将样本序列制成随机游走（random walk）； $g$ : 假定密度/跳跃分布(jumping distribution)

2. 迭代 $t$

通过密度函数 $g(x'|x_t)$ 生成备选样本值 $x'$
计算接受备选样本值的比率: $\alpha=f(x')/f(x_t)$ ；由于 $f(x)$ 又与 $P(x)$ 成比例，因此 $\alpha=f(x')/f(x_t)=P(x')/P(x_t)$ 。
若 $\alpha \geq 1$ ， $x'$ 比 $x_t$ 更接近 $P(x)$ ，令 $x_{t+1}=x'$ ；否则，以 $\alpha$ 的概率接受备选样本值 $x'$ ；若备选样本值被拒绝，令 $x_{t+1}=x_t$ 。

1.5 缺陷

样本自相关：相邻的样本会相互相关，虽然可以通过每 $n$ 步取样的方式来减少相关性，但这样的后果就是很难让样本近似于目标分布 $P(x)$ 。

自相关性可通过增加步调长度（jumping width，与jumping function $g(x|y)$ 的方差有关）来控制，但同时也增加了拒绝备选样本的几率 $\alpha$ 。
过大或过小的jumping size会导致slow-mixing Markov Chain, 即高度自相关的一组样本，以至于我们需要得到非常大的样本量 $n$ 才能得到目标分布 $P(x)$ 。

初始值的选择：尽管Markov Chain最后都会收敛到目标分布，然而初始值的选择直接影响到运算时间，尤其是把初始值选在在了“低密度”区域。因此，选择初值时，最好加入一个“burn-in period”（预烧期，预选期）。

1.6 优势

抗“高维魔咒”（curse of dimensionality）：维度增加，对于rejection sampling方法来说，拒绝的概率就是呈指数增长。而MCMC则成为了解决这种问题的唯一方法
多元分布中，为了避免多元初始值以及 $g(x|y)$ 选择不当而导致的问题，Gibbs sampling是另外一个更适合解决多元分布问题的MCMC 方法。Gibbs sampling从多元分布的各个维度中分别选择初始值，然后这些变量分别同时取样。

1.7 衍生算法

Metropolis-Hastings算法的目的是根据目标分布 $P(x)$ 生成一个状态集。算法运用了渐进于平稳分布 $\pi(x)$ 的Markov过程使得 $\pi(x)=P(x)$ 。

Markov过程由它的状态转换概率所决定的。 $P(x'|x)$ 为从任意状态 $x$ 转化到另外一个任意状态 $x'$ 的转换概率。当它符合以下两点是，它有唯一的平稳分布 $\pi(x)$ :

平稳分布 $\pi(x)$ 的存在性：充分但非必要条件为细致平衡（detailed balance）：每一个转换 $x\rightarrow x'$ 都是可逆的，即对任意连个状态 $x$ 和 $x'$ ，在状态 $x$ 时 $x \rightarrow x'$ 的概率等于在状态 $x'$ 时 $x' \rightarrow x$ 的概率，即 $\pi(x) P(x'|x)=\pi(x')P(x|x')$
平稳分布 $\pi(x)$ 的唯一性：由Markov过程的遍历性可得，即每一个转换 $x\rightarrow x'$ 需满足：

非周期性的（aperiodic）：系统不会再某一个固定时间区间内回到原状态
正递归的（positive recurrent）：回到同一个状态的步数的期望是有限值

Metropolis-Hastings算法需要通过构建转换概率来构建一个Markov过程且需要满足以上两点来使得Markov过程的平稳分布 $\pi(x)$ 就是目标分布 $P(x)$ 。算法衍生从满足细致平衡条件开始：

$P(x'|x)P(x)=P(x|x')P(x')$ ,

即 $\dfrac{P(x'|x)}{P(x|x')}=\dfrac{P(x')}{P(x)}$ (1)

这样就把转换的细致平衡转化为两步：提议和接受-拒绝。

提议分布 $g(x'|x)$ 为给定状态 $x$ 时提出状态 $x'$ 的条件概率，接受分布 $A(x'|x)$ 为接受所提出的状态 $x'$ 的概率。转换概率就能写成它们的乘积：

$P(x'|x)=g(x'|x)A(x'|x)$ (2)

由(1)和(2)就有：

$\dfrac{A(x'|x)}{A(x|x')}=\dfrac{P(x')}{P(x)} \dfrac{g(x|x')}{g(x'|x)}$

之后就要选一个满足上面条件的接受函数 $A(x'|x)$ 了。一般我们选Metropolis选项，即：

$A(x'|x)=min \left(1, \dfrac{P(x')}{P(x)} \dfrac{g(x|x')}{g(x'|x)} \right)$

满足了Metropolis算法的基本步骤，即大于1的时候总是接受，小于1的时候按比例接受和拒绝。

则Metropolis-Hastings算法可以概述为：

初始化：随机选取初始状态 $x$
根据 $g(x'|x)$ 随机选取状态 $x'$
根据 $A(x'|x)$ 选择接受或拒绝状态 $x'$ ；若接受，系统转换到状态 $x'$ ，否则停留在状态 $x$
返回步骤2，直到生成了 $T$ 个状态
保存状态 $x$ ，返回步骤2

原则上说，被保存的状态是从分布 $P(x)$ 中生成的，因为步骤4保证了这些状态是非相关的。取值 $T$ 需根据不同的特性选取，如提出的分布以及，它必须是Markov过程自相关时间的阶。

关于 $g(x'|x)$ 的选取需具体问题具体分析和调整。

先验：

 
epsilon<-pars[3]    #残差
 
 
 
prior_a<-dnorm(a,0,100,log=TRUE)     ##所有的非信息性先验
 
prior_b<-dnorm(b,0,100,log=TRUE)     ## 参数.
 
prior_epsilon<-dgamma(epsilon,1,1/100,log=TRUE)

1）现在让我们模拟一些数据以进行运行测试：

x<-runif(30,5,15)
 
y<-x+rnorm(30,0,5) ##斜率=1, 截距=0, epsilon=5

2）Metro Hastings 完成所有工作。

MH(li_func=li_reg,pars=c(0,1,1),

视频

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

探索见解 ➜

去bilibili观看 ➜

探索更多视频 ➜

3）您可以使用plotMH（）查看所有模型参数的后验

plot(mcmc)

绘制所有参数之间的相关性。

4）输出后验置信区间。

BCI
 
#              0.025    0.975
 
# a       -5.3345970 6.841016
 
# b        0.4216079 1.690075
 
# epsilon  3.8863393 6.660037

接下来，我想提供一种直观的方法来可视化此算法运行的情况。

主要思想是从分布中抽取样本。积分很重要，贝叶斯定理本身：

P（θ| D）= P（D |θ）P（θ）/ P（D）

其中P（D）是观察数据的无条件概率。由于这不依赖于推断的模型（θ）参数，因此P（D）是归一化常数。

因此，我们有一个非归一化的概率密度函数，我们希望通过随机抽样来估计。对于复杂的模型而言，随机抽样本身的过程通常很困难，因此，我们使用马尔可夫链来探索分布。我们需要一个链，如果运行时间足够长，它将作为目标分布的随机样本整体。我们构建的马尔可夫链的这种特性称为 遍历性。Metropolis-Hastings算法是构建这种链的一种方法。

步骤

在参数空间k_X中选择一些起点
选择一个候选点k_Y〜N（k_X，σ）。这通常称为提议分布。
移至候选点的概率为：min（π（k_Y）/π（K_X），1）
重复。

最受欢迎的见解

1.使用R语言进行METROPLIS-IN-GIBBS采样和MCMC运行

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

3.R语言实现MCMC中的Metropolis–Hastings算法与吉布斯采样

4.R语言BUGS JAGS贝叶斯分析马尔科夫链蒙特卡洛方法（MCMC）采样

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

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

7.R语言用Rcpp加速Metropolis-Hastings抽样估计贝叶斯逻辑回归模型的参数

8.R语言使用Metropolis- Hasting抽样算法进行逻辑回归

9.R语言中基于混合数据抽样(MIDAS)回归的HAR-RV模型预测GDP增长

以下代码通过简单的正态目标分布演示了此过程。

 
###     Metropolis-Hastings 可视化                #######
 
 
 
k_X = seed; ##将k_X设置为种子位置
 
 
 
 
for(i in 1:iter)
 
{
 
track<-c(track,k_X)    ## 链
 
k_Y = rnorm(1,k_X,prop_sd) ##候选点
 
 
 
## -- 绘制链的核密度估计
 
 
 
 
 
lines(density(track,adjust=1.5),col='red',lwd=2)
 
 
 
## -- 绘制链
 
 
plot(track,1:i,xlim=plot_range,main='',type='l',ylab='Trace')
 
 
 
 
## -- 绘制目标分布和提议分布 
 
 
 
curve(dnorm(x,k_X,prop_sd),col='black',add=TRUE)
 
abline(v=k_X,lwd=2)
 
 
 
 
## 接受概率为a_X_Y 
 
if (log(runif(1))<=a_X_Y)
 
 
 
points(k_Y,0,pch=19,col='green',cex=2)
 
 
 
## 调整提议
 
if(i>100)
 
prop_sd=sd(track[floor(i/2):i])