共计 2685 个字符，预计需要花费 7 分钟才能阅读完成。
原文链接：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.3678794target(-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 理解以下内容：
- 不同的起始值如何影响 MCMC 计划？
- 较大 / 较小的提案标准差有什么影响？
- 尝试将指标函数更改为
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 增长
