R语言实现MCMC中的Metropolis–Hastings算法与吉布斯采样
原文连接：http://tecdat.cn/?p=3772
建立测试数据
第一步，咱们建立一些测试数据，用来拟合咱们的模型。咱们假设预测变量和因变量之间存在线性关系，因此咱们用线性模型并添加一些噪音。python
trueA <- 5 trueB <- 0 trueSd <- 10 sampleSize <- 31 # 建立独立的x值 x <- (-(sampleSize-1)/2):((sampleSize-1)/2) # 根据ax + b + N（0，sd）建立因变量 y <- trueA * x + trueB + rnorm(n=sampleSize,mean=0,sd=trueSd) plot(x,y, main="Test Data")
图算法
定义统计模型
下一步是指定统计模型。咱们已经知道数据是用x和y之间的线性关系y = a * x + b和带有标准差sd的正态偏差模型N（0，sd）建立的，因此让咱们使用相同的模型进行拟合，看看若是咱们能够检索咱们的原始参数值。编程
从模型中导出似然函数
为了估计贝叶斯分析中的参数，咱们须要导出咱们想要拟合的模型的似然函数。似然函数是咱们指望观察到的数据以咱们所看到的模型的参数为条件发生的几率（密度）。所以，鉴于咱们的线性模型y = b + a*x + N(0,sd)将参数（a, b, sd）做为输入，咱们必须返回在这个模型下得到上述测试数据的几率（这听起来比较复杂，正如你在代码中看到的，咱们只是计算预测值y = b + a*x与观察到的y之间的差别，而后咱们必须查找这种误差发生的几率密度（使用dnorm）。函数
likelihood <- function(param){
a = param\[1\]
b = param\[2\]
sd = param\[3\]
pred = a*x + b
sumll = sum(singlelikelihoods)
(sumll)
}
slopevalues <- function(x){return(likelihood(c(x, trueB, trueSd)))}
斜率参数对数似然曲线post
做为说明，代码的最后几行绘制了斜率参数a的一系列参数值的似然函数。学习
为何咱们使用对数
您注意到结果是似然函数中几率的对数，这也是我对全部数据点的几率求和的缘由（乘积的对数等于对数之和）。咱们为何要作这个？ 由于不少小几率乘以的可能性很快就会变得很是小（好比10 ^ -34）。在某些阶段，计算机程序存在数字四舍五入的问题。 测试
定义先验
第二步，与贝叶斯统计中同样，咱们必须为每一个参数指定先验分布。为了方便起见，我对全部三个参数使用了均匀分布和正态分布。 _无信息先验一般是1 / sigma（若是你想了解缘由，请看这里）。 _优化
#先验分布
prior <- function(param){
a = param\[1\]
b = param\[2\]
sd = param\[3\]
aprior = (a, min=0, max=10, log = T)
bprior = dnorm(b, sd = 5, log = T)
}
后验
先验和似然性的乘积是MCMC将要处理的实际数量。这个函数被称为后验 。一样，在这里咱们使用加总，由于咱们使用对数。spa
posterior <- function(param){
return ( (param) + prior(param))
}
MCMC
接下来是Metropolis-Hastings算法。该算法最多见的应用之一（如本例所示）是从贝叶斯统计中的后验密度中提取样本。然而，原则上，该算法可用于从任何可积函数中进行采样。所以，该算法的目的是在参数空间中跳转，可是以某种方式使得在某一点上的几率与咱们采样的函数成比例（这一般称为目标函数）。在咱们的例子中，这是上面定义的后验。code
这是经过
- 从随机参数值开始
- 根据称为提议函数的某个几率密度，选择接近旧值的新参数值
- 以几率p（新）/ p（旧）跳到这个新点，其中p是目标函数，p> 1表示跳跃
当咱们运行这个算法时，它访问的参数的分布会收敛到目标分布p。那么，让咱们在R中获得 ：
########Metropolis算法# ################
proposalfunction <- function(param){
return(rnorm(3,mean = param, sd= c(0.1,0.5,0.3)))
}
run\_metropolis\_MCMC <- function(startvalue, iterations){
for (i in 1:iterations){
if (runif(1) < probab){
chain\[i+1,\] = proposal
}else{
chain\[i+1,\] = chain\[i,\]
}
}
return(chain)
}
chain = run\_metropolis\_MCMC(startvalue, 10000)
burnIn = 5000
acceptance = 1-mean(duplicated(chain\[-(1:burnIn),\]))
使用后验的对数可能在开始时有点混乱，特别是当您查看计算接受几率的行时（probab = exp（后验分布（提议分布） - 后验（链[i，]）） ）。要理解咱们为何这样作，请注意p1 / p2 = exp [log（p1）-log（p2）]。
算法的第一步可能受初始值的误差，所以一般被丢弃用于进一步分析 。要看的一个有趣的输出是接受率： 接受标准拒绝提议的频率是多少？接受率能够受提议函数的影响：一般，提议越接近，接受率越大。然而，很是高的接受率一般是无益的：这意味着算法“停留”在同一点 。能够证实，20％到30％的接受率对于典型应用来讲是最佳的 。
###概要####################### par(mfrow = c(2,3)) hist( \[-(1:burnIn),1\],nclass=30, , main="Posterior of a", xlab="True value = red line" ) abline(v = mean(chain\[-(1:burnIn),1\])) #进行比较： summary(lm(y~x))
生成的图应该相似于下图。您能够看到咱们或多或少地检索了用于建立数据的原始参数，而且您还看到咱们得到了一个围绕最高后验值的特定区域，这些区域也有一些数据支持，这是贝叶斯至关于置信度的区间。
所获得的图应该看起来像下面的图。你看到咱们检索到了或多或少用于建立数据的原始参数，你还看到咱们在最高后验值周围获得了必定的区域，这些后验值也有一些数据，这至关于贝叶斯的置信区间。
图：上排显示斜率（a）的后验估计，截距（b）和偏差的标准误差（sd）。下一行显示马尔可夫链参数值。
还有问题吗？请在下面留言！
最受欢迎的看法
4.R语言中的block Gibbs吉布斯采样贝叶斯多元线性回归