原文链接：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     #这只是一个起始值，我设置为3for(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增长