乐趣区

关于算法:R语言实现MCMC中的Metropolis–Hastings算法与吉布斯采样

jiezi

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

创立测试数据

第一步，咱们创立一些测试数据，用来拟合咱们的模型。咱们假如预测变量和因变量之间存在线性关系，所以咱们用线性模型并增加一些乐音。

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)))}

斜率参数对数似然曲线

作为阐明，代码的最初几行绘制了斜率参数 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 将要解决的理论数量。这个函数被称为后验。同样，在这里咱们应用加总，因为咱们应用对数。

posterior <- function(param){return ( (param) + prior(param))

}

MCMC

接下来是 Metropolis-Hastings 算法。该算法最常见的利用之一（如本例所示）是从贝叶斯统计中的后验密度中提取样本。然而，原则上，该算法可用于从任何可积函数中进行采样。因而，该算法的目标是在参数空间中跳转，然而以某种形式使得在某一点上的概率与咱们采样的函数成比例（这通常称为指标函数）。在咱们的例子中，这是下面定义的后验。

这是通过

  1. 从随机参数值开始
  2. 依据称为提议函数的某个概率密度，抉择靠近旧值的新参数值
  3. 以概率 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）。下一行显示马尔可夫链参数值。

还有问题吗？请在上面留言！

最受欢迎的见解

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 模型实现

退出移动版