关于数据挖掘:R语言贝叶斯推断与MCMC实现MetropolisHastings-采样算法示例附代码数据

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

原文出处：拓端数据部落公众号

 
最近咱们被客户要求撰写对于贝叶斯推断的钻研报告，包含一些图形和统计输入。

示例1：应用MCMC的指数分布采样

任何MCMC计划的指标都是从“指标”散布产生样本。在这种状况下，咱们将应用平均值为1的指数分布作为咱们的指标散布。所以咱们从定义指标密度开始：

target = function(x){
  if(x<0){
    return(0)}
  else {
    return( exp(-x))
  }
}

定义了函数之后，咱们当初能够用它来计算几个值（只是为了阐明函数的概念）：

target(1)

[1] 0.3678794

target(-1)

[1] 0

接下来，咱们将布局一个Metropolis-Hastings计划，从与指标成比例的散布中进行抽样

x[1] = 3     #这只是一个起始值，我设置为3
for(i in 2:1000){
  A = target(proposedx)/target(currentx) 
  if(runif(1)<A){
    x[i] = proposedx       # 承受概率min（1，a）
  } else {
    x[i] = currentx        #否则“回绝”口头，放弃原样
  }

留神，x是马尔可夫链的实现。咱们能够画几个x的图：

咱们能够将其封装在一个mcmc函数中，以使代码更整洁，这样更改起始值和提议散布更容易

  for(i in 2:niter){
    currentx = x[i-1]
    proposedx = rnorm(1,mean=currentx,sd=proposalsd) 
    A = target(proposedx)/target(currentx)
    if(runif(1)<A){
      x[i] = proposedx       # 承受概率min（1，a）
    } else {
      x[i] = currentx        # 否则“回绝”口头，放弃原样
    }

当初咱们将运行MCMC计划3次，看看后果有多类似：

z1=MCMC(1000,3,1)
z2=MCMC(1000,3,1)
z3=MCMC(1000,3,1)

plot(z1,type="l")

par(mfcol=c(3,1)) #通知R将3个图形放在一个页面上

hist(z1,breaks=seq(0,maxz,length=20))

练习

应用函数easyMCMC理解以下内容：

1. 不同的起始值如何影响MCMC计划？
2. 较大/较小的提案标准差有什么影响？
3. 尝试将指标函数更改为

target = function(x){
  
  return((x>0 & x <1) + (x>2 & x<3))
}

这个指标看起来像什么？如果提议sd太小怎么办？（例如，尝试1和0.1）

例2：预计等位基因频率

在对双等位基因座的基因型（如具备AA和AA等位基因的基因座）进行建模时，一个规范的假如是群体是“随机”的。这意味着如果p是等位基因AA的频率，那么基因型和将别离具备频率和。

p一个简略的先验是假如它在[0,1]上是平均的。 假如咱们抽样n个个体，察看基因型、基因型和基因型。

上面的R代码给出了一个简短的MCMC例程，能够从p的后验散布中进行采样。请尝试遍历该代码，看看它是如何工作的。

prior = function(p){
  if((p<0) || (p>1)){  # || 这里意思是“或”
    return(0)}
  else{
    return(1)}
}

likelihood = function(p, nAA, nAa, naa){
  return(p^(2*nAA) * (2*p*(1-p))^nAa * (1-p)^(2*naa))
}

psampler = function(nAA, nAa, naa){

  for(i in 2:niter){

    if(runif(1)<A){
      p[i] = newp       # 承受概率min（1，a）
    } else {
      p[i] = currentp        # 否则“回绝”口头，放弃原样
    }

对运行此样本。

当初用一些R代码来比拟后验样本和实践后验样本（在这种状况下能够通过剖析取得；因为咱们察看到121个As和79个as，在200个样本中，p的后验样本是β（121+1,79+1）。

hist(z,prob=T)
lines(x,dbeta(x,122, 80))  # 在直方图上叠加β密度

您也可能心愿将前5000 z的值抛弃为“burnin”（预烧期）。这里有一种办法，在R中仅抉择最初5000 z

hist(z[5001:10000])

练习

钻研终点和提案标准偏差如何影响算法的收敛性。

例3：预计等位基因频率和近交系数

一个简单一点的代替办法是假如人们有一种偏向，即人们会与比“随机”关系更亲密的其他人近交（例如，在天文构造上的人口中可能会产生这种状况）。一个简略的办法是引入一个额定的参数，即“近亲繁殖系数”f，并假如和基因型有频率

和。

在大多数状况下，将f作为种群特色来看待是很天然的，因而假如f在各个位点上是恒定的。

请留神，f和p都被束缚为介于0和1之间（包含0和1）。对于这两个参数中的每一个，一个简略的先验是假如它们在[0,1]上是独立的。 假如咱们抽样n个个体，察看基因型、基因型和基因型。

练习：

• 编写一个短的MCMC程序，从f和p的联结散布中取样。

sampler = function(){

  f[1] = fstartval
  p[1] = pstartval
  for(i in 2:niter){
    currentf = f[i-1]
    currentp = p[i-1]
    newf = currentf + 
    newp = currentp + 
  
  }
  return(list(f=f,p=p)) # 返回一个蕴含两个名为f和p的元素的“list”
}

• 应用此样本取得f和p的点估计（例如，应用后验平均数）和f和p的区间预计（例如，90%后验置信区间），数据：

附录：GIBBS采样

您也能够用Gibbs采样器解决这个问题 

为此，您将想要应用以下“潜在变量”示意模型：

将zi相加失去与上述雷同的模型:

最受欢迎的见解

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增长