乐趣区

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

jiezi

原文链接：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 增长

退出移动版