原文链接：http://tecdat.cn/?p=19889

如果您可以写出模型的似然函数，则 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算法是构建这种链的一种方法。

步骤：