0
点赞
收藏
分享

微信扫一扫

拓端数据tecdat|R语言贝叶斯推断与MCMC:实现Metropolis-Hastings 采样算法示例

践行数据分析 2022-11-14 阅读 147

 

示例1:使用MCMC的指数分布采样

任何MCMC方案的目标都是从“目标”分布产生样本。在这种情况下,我们将使用平均值为1的指数分布作为我们的目标分布。所以我们从定义目标密度开始:

  1.  target = function(x){
  2.  if(x<0){
  3.  return(0)}
  4.  else {
  5.  return( exp(-x))
  6.  }
  7.  }

定义了函数之后,我们现在可以用它来计算几个值(只是为了说明函数的概念):

target(1)

[1] 0.3678794

target(-1)

[1] 0

接下来,我们将规划一个Metropolis-Hastings方案,从与目标成比例的分布中进行抽样

  1.  x[1] = 3 #这只是一个起始值,我设置为3
  2.  for(i in 2:1000){
  3.  A = target(proposedx)/target(currentx)
  4.  if(runif(1)<A){
  5.  x[i] = proposedx # 接受概率min(1,a)
  6.  } else {
  7.  x[i] = currentx #否则“拒绝”行动,保持原样
  8.  }

注意,x是马尔可夫链的实现。我们可以画几个x的图:

拓端数据tecdat|R语言贝叶斯推断与MCMC:实现Metropolis-Hastings 采样算法示例_数据

 

 

拓端数据tecdat|R语言贝叶斯推断与MCMC:实现Metropolis-Hastings 采样算法示例_r语言_02

 

 

我们可以将其封装在一个mcmc函数中,以使代码更整洁,这样更改起始值和提议分布更容易

  1.   
  2.  for(i in 2:niter){
  3.  currentx = x[i-1]
  4.  proposedx = rnorm(1,mean=currentx,sd=proposalsd)
  5.  A = target(proposedx)/target(currentx)
  6.  if(runif(1)<A){
  7.  x[i] = proposedx # 接受概率min(1,a)
  8.  } else {
  9.  x[i] = currentx # 否则“拒绝”行动,保持原样
  10.  }

现在我们将运行MCMC方案3次,看看结果有多相似:

  1.  z1=MCMC(1000,3,1)
  2.  z2=MCMC(1000,3,1)
  3.  z3=MCMC(1000,3,1)
  4.   
  5.  plot(z1,type="l")

拓端数据tecdat|R语言贝叶斯推断与MCMC:实现Metropolis-Hastings 采样算法示例_r语言_03

 

 

  1.  par(mfcol=c(3,1)) #告诉R将3个图形放在一个页面上
  2.   
  3.  hist(z1,breaks=seq(0,maxz,length=20))

拓端数据tecdat|R语言贝叶斯推断与MCMC:实现Metropolis-Hastings 采样算法示例_数据_04

 

 

练习

使用函数easyMCMC了解以下内容:

  1. 不同的起始值如何影响MCMC方案?
  2. 较大/较小的提案标准差有什么影响?
  3. 尝试将目标函数更改为
  4.  target = function(x){
  5.   
  6.  return((x>0 & x <1) + (x>2 & x<3))
  7.  }

这个目标看起来像什么?如果提议sd太小怎么办?(例如,尝试1和0.1)

例2:估计等位基因频率

在对双等位基因座的基因型(如具有AA和AA等位基因的基因座)进行建模时,一个标准的假设是群体是“随机”的。这意味着如果p是等位基因AA的频率,那么基因型

拓端数据tecdat|R语言贝叶斯推断与MCMC:实现Metropolis-Hastings 采样算法示例_sed_05

拓端数据tecdat|R语言贝叶斯推断与MCMC:实现Metropolis-Hastings 采样算法示例_数据_06

将分别具有频率

拓端数据tecdat|R语言贝叶斯推断与MCMC:实现Metropolis-Hastings 采样算法示例_数据_07


拓端数据tecdat|R语言贝叶斯推断与MCMC:实现Metropolis-Hastings 采样算法示例_r语言_08

。p一个简单的先验是假设它在[0,1]上是均匀的。 假设我们抽样n个个体,观察

拓端数据tecdat|R语言贝叶斯推断与MCMC:实现Metropolis-Hastings 采样算法示例_r语言_09

基因型

拓端数据tecdat|R语言贝叶斯推断与MCMC:实现Metropolis-Hastings 采样算法示例_sed_10


拓端数据tecdat|R语言贝叶斯推断与MCMC:实现Metropolis-Hastings 采样算法示例_sed_11

基因型

拓端数据tecdat|R语言贝叶斯推断与MCMC:实现Metropolis-Hastings 采样算法示例_r语言_12

拓端数据tecdat|R语言贝叶斯推断与MCMC:实现Metropolis-Hastings 采样算法示例_sed_13基因型

拓端数据tecdat|R语言贝叶斯推断与MCMC:实现Metropolis-Hastings 采样算法示例_r语言_14


下面的R代码给出了一个简短的MCMC例程,可以从p的后验分布中进行采样。请尝试遍历该代码,看看它是如何工作的。

  1.  prior = function(p){
  2.  if((p<0) || (p>1)){ # || 这里意思是“或”
  3.  return(0)}
  4.  else{
  5.  return(1)}
  6.  }
  7.   
  8.  likelihood = function(p, nAA, nAa, naa){
  9.  return(p^(2*nAA) * (2*p*(1-p))^nAa * (1-p)^(2*naa))
  10.  }
  11.   
  12.  psampler = function(nAA, nAa, naa){
  13.   
  14.  for(i in 2:niter){
  15.   
  16.  if(runif(1)<A){
  17.  p[i] = newp # 接受概率min(1,a)
  18.  } else {
  19.  p[i] = currentp # 否则“拒绝”行动,保持原样
  20.  }


拓端数据tecdat|R语言贝叶斯推断与MCMC:实现Metropolis-Hastings 采样算法示例_r语言_15

运行此样本。

现在用一些R代码来比较后验样本和理论后验样本(在这种情况下可以通过分析获得;因为我们观察到121个As和79个as,在200个样本中,p的后验样本是β(121+1,79+1)。

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

拓端数据tecdat|R语言贝叶斯推断与MCMC:实现Metropolis-Hastings 采样算法示例_数据_16

您也可能希望将前5000 z的值丢弃为“burnin”(预烧期)。这里有一种方法,在R中仅选择最后5000 z

hist(z[5001:10000])

拓端数据tecdat|R语言贝叶斯推断与MCMC:实现Metropolis-Hastings 采样算法示例_r语言_17

 

 

练习

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

例3:估计等位基因频率和近交系数

一个复杂一点的替代方法是假设人们有一种倾向,即人们会与比“随机”关系更密切的其他人近交(例如,在地理结构上的人口中可能会发生这种情况)。一个简单的方法是引入一个额外的参数,即“近亲繁殖系数”f,并假设

拓端数据tecdat|R语言贝叶斯推断与MCMC:实现Metropolis-Hastings 采样算法示例_数据_18


拓端数据tecdat|R语言贝叶斯推断与MCMC:实现Metropolis-Hastings 采样算法示例_数据_19

基因型有频率

拓端数据tecdat|R语言贝叶斯推断与MCMC:实现Metropolis-Hastings 采样算法示例_数据_20


拓端数据tecdat|R语言贝叶斯推断与MCMC:实现Metropolis-Hastings 采样算法示例_sed_21


在大多数情况下,将f作为种群特征来对待是很自然的,因此假设f在各个位点上是恒定的。

请注意,f和p都被约束为介于0和1之间(包括0和1)。对于这两个参数中的每一个,一个简单的先验是假设它们在[0,1]上是独立的。 假设我们抽样n个个体,观察

拓端数据tecdat|R语言贝叶斯推断与MCMC:实现Metropolis-Hastings 采样算法示例_r语言_22

基因型

拓端数据tecdat|R语言贝叶斯推断与MCMC:实现Metropolis-Hastings 采样算法示例_sed_23


拓端数据tecdat|R语言贝叶斯推断与MCMC:实现Metropolis-Hastings 采样算法示例_r语言_24

基因型

拓端数据tecdat|R语言贝叶斯推断与MCMC:实现Metropolis-Hastings 采样算法示例_r语言_25

拓端数据tecdat|R语言贝叶斯推断与MCMC:实现Metropolis-Hastings 采样算法示例_r语言_26基因型

拓端数据tecdat|R语言贝叶斯推断与MCMC:实现Metropolis-Hastings 采样算法示例_r语言_27


练习:

  • 编写一个短的MCMC程序,从f和p的联合分布中取样。
  1.  sampler = function(){
  2.   
  3.  f[1] = fstartval
  4.  p[1] = pstartval
  5.  for(i in 2:niter){
  6.  currentf = f[i-1]
  7.  currentp = p[i-1]
  8.  newf = currentf +
  9.  newp = currentp +
  10.   
  11.  }
  12.  return(list(f=f,p=p)) # 返回一个包含两个名为f和p的元素的“list”
  13.  }
  • 使用此样本获得f和p的点估计(例如,使用后验平均数)和f和p的区间估计(例如,90%后验置信区间),数据:

 

附录:GIBBS采样

您也可以用Gibbs采样器解决这个问题 

为此,您将想要使用以下“潜在变量”表示模型:

拓端数据tecdat|R语言贝叶斯推断与MCMC:实现Metropolis-Hastings 采样算法示例_数据_28

将zi相加得到与上述相同的模型:

 

拓端数据tecdat|R语言贝叶斯推断与MCMC:实现Metropolis-Hastings 采样算法示例_r语言_29

拓端数据tecdat|R语言贝叶斯推断与MCMC:实现Metropolis-Hastings 采样算法示例_sed_30


举报

相关推荐

0 条评论