示例 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 相加失去与上述雷同的模型:

 ​