R语言贝叶斯推断与MCMC:实现Metropolis-Hastings 采样算法示例

简介: R语言贝叶斯推断与MCMC:实现Metropolis-Hastings 采样算法示例

示例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语言Stan贝叶斯回归置信区间后验分布可视化模型检验|附数据代码
R语言Stan贝叶斯回归置信区间后验分布可视化模型检验|附数据代码
|
1月前
|
数据可视化 定位技术
R语言贝叶斯INLA空间自相关、混合效应、季节空间模型、SPDE、时空分析野生动物数据可视化
R语言贝叶斯INLA空间自相关、混合效应、季节空间模型、SPDE、时空分析野生动物数据可视化
|
1月前
|
存储 机器学习/深度学习 算法
R语言贝叶斯Metropolis-Hastings采样 MCMC算法理解和应用可视化案例
R语言贝叶斯Metropolis-Hastings采样 MCMC算法理解和应用可视化案例
|
1月前
|
数据采集 机器学习/深度学习 数据可视化
R语言贝叶斯模型预测电影评分数据可视化分析
R语言贝叶斯模型预测电影评分数据可视化分析
|
3天前
|
算法 JavaScript 决策智能
基于禁忌搜索算法的TSP路径规划matlab仿真
**摘要:** 使用禁忌搜索算法解决旅行商问题(TSP),在MATLAB2022a中实现路径规划,显示优化曲线与路线图。TSP寻找最短城市访问路径,算法通过避免局部最优,利用禁忌列表不断调整顺序。关键步骤包括初始路径选择、邻域搜索、解评估、选择及禁忌列表更新。过程示意图展示搜索效果。
|
3天前
|
机器学习/深度学习 算法
基于BP神经网络和小波变换特征提取的烟草香型分类算法matlab仿真,分为浓香型,清香型和中间香型
```markdown 探索烟草香型分类:使用Matlab2022a中的BP神经网络结合小波变换。小波分析揭示香气成分的局部特征,降低维度,PCA等用于特征选择。BP网络随后处理这些特征,以区分浓香、清香和中间香型。 ```
|
4天前
|
机器学习/深度学习 算法
m基于PSO-GRU粒子群优化长门控循环单元网络的电力负荷数据预测算法matlab仿真
摘要: 在MATLAB 2022a中,对比了电力负荷预测算法优化前后的效果。优化前为&quot;Ttttttt111222&quot;,优化后为&quot;Tttttttt333444&quot;,明显改进体现为&quot;Tttttttttt5555&quot;。该算法结合了粒子群优化(PSO)和长门控循环单元(GRU)网络,利用PSO优化GRU的超参数,提升预测准确性和稳定性。PSO模仿鸟群行为寻找最优解,而GRU通过更新门和重置门处理长期依赖问题。核心MATLAB程序展示了训练和预测过程,包括使用&#39;adam&#39;优化器和超参数调整,最终评估并保存预测结果。
13 0
|
5天前
|
算法 安全
基于龙格库塔算法的SIR病毒扩散预测matlab仿真
该程序使用龙格库塔算法实现SIR模型预测病毒扩散,输出易感、感染和康复人群曲线。在MATLAB2022a中运行显示预测结果。核心代码设置时间区间、参数,并定义微分方程组,通过Runge-Kutta方法求解。SIR模型描述三类人群动态变化,常微分方程组刻画相互转化。模型用于预测疫情趋势,支持公共卫生决策,但也存在局限性,如忽略空间结构和人口异质性。
|
5天前
|
机器学习/深度学习 监控 算法
基于yolov2深度学习网络的昆虫检测算法matlab仿真,并输出昆虫数量和大小判决
YOLOv2算法应用于昆虫检测,提供实时高效的方法识别和定位图像中的昆虫,提升检测精度。核心是统一检测网络,预测边界框和类别概率。通过预测框尺寸估算昆虫大小,适用于农业监控、生态研究等领域。在matlab2022A上运行,经过关键升级,如采用更优网络结构和损失函数,保证速度与精度。持续优化可增强对不同昆虫的检测能力。![image.png](https://ucc.alicdn.com/pic/developer-ecology/3tnl7rfrqv6tw_e760ff6682a3420cb4e24d1e48b10a2e.png)
|
6天前
|
算法 调度 决策智能
基于自适应遗传算法的车间调度matlab仿真,可以任意调整工件数和机器数,输出甘特图
这是一个使用MATLAB2022a实现的自适应遗传算法解决车间调度问题的程序,能调整工件数和机器数,输出甘特图和适应度收敛曲线。程序通过编码初始化、适应度函数、遗传操作(选择、交叉、变异)及自适应机制进行优化,目标如最小化完工时间。算法在迭代过程中动态调整参数,以提升搜索效率和全局优化。

热门文章

最新文章