如果您可以写出模型的似然函数，则 Metropolis-Hastings算法可以负责其余部分（即MCMC ）。我写了r代码来简化对任意模型的后验分布的估计。具体如下：
1）定义模型（即概率先验）。在此示例中，让我们构建一个简单的线性回归模型（对数）。
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) )
先验：
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)
现在让我们模拟一些数据以进行运行测试：
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),
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）
- 重复。
以下代码通过简单的正态目标分布演示了此过程。
### 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])
该算法实现中的一个普遍问题是σ的选择。当σ接近目标分布的标准偏差时，将发生有效混合（链收敛到目标分布）。当我们不知道这个值时。我们可以允许σ根据到目前为止的链历史记录进行调整。在上面的示例中，将σ更新为链中某些先验点的标准偏差值。
输出为多页pdf，可以滚动浏览。
顶部显示了目标分布（蓝色虚线）和通过MCMC样本对目标进行的核平滑估计。第二面板显示了链的轨迹，底部显示了算法本身的步骤。
注意：请注意，前100次左右的迭代是目标分布的较差表示。在实践中，我们将“预烧”该链的前n个迭代-通常是前100-1000个。
