原文链接:http://tecdat.cn/?p=24721  

本文,我通过两个种群生态学家可能感兴趣的例子来阐明应用“JAGS”来模仿数据:首先是线性回归,其次是预计动物存活率(公式化为状态空间模型)。

最近,我始终在致力模仿来自简单分层模型的数据。我当初正在应用 JAGS

模仿数据 JAGS 很不便,因为你能够应用(简直)雷同的代码进行模仿和推理,并且你能够在雷同的环境(即JAGS)中进行模仿钻研(偏差、精度、区间笼罩 )。

线性回归示例

咱们首先加载本教程所需的包:

library(R2jags)

而后间接切入正题,让咱们从线性回归模型生成数据。应用一个 data 块,并将参数作为数据传递。

data{# 似然函数:for (i in 1:N){y\[i\] ~  # tau是精度(1/方差)。}

这里, alpha 和 beta 是截距和斜率、 tau 方差的精度或倒数、 y 因变量和 x 解释变量。

咱们为用作数据的模型参数抉择一些值:

# 模仿的参数 N  # 样本x <- 1:N # 预测因子alpha # 截距beta  # 斜率sigma# 残差sd 1/(sigma*sigma) # 精度# 在模仿步骤中,参数被当作数据处理

当初运行 JAGS; 请留神,咱们监控因变量而不是参数,就像咱们在进行规范推理时所做的那样:

# 运行后果out

输入有点乱,须要适当格式化:

# 从新格式化输入mcmc(out)

dim

dat

当初让咱们将咱们用来模仿的模型拟合到咱们刚刚生成的数据中。不再赘述,假如读者相熟 JAGS 线性回归。

# 用BUGS语言指定模型model <-     for (i in 1:N){y\[i\] ~ dnorm(mu\[i\], tau) # tau是精度(1/方差)alpha  截距beta # 斜率sigma  # 标准差# 数据dta <- list(y = dt, N = length(at), x = x)# 初始值inits # MCMC设置ni <- 10000# 从R中调用JAGSjags()

让咱们看看后果并与咱们用来模仿数据的参数进行比拟(见上文):

# 总结后验print(res)

查看收敛:

# 追踪图plot(res)

绘制回归参数和残差标准差的后验散布:

# 后验散布plot(res)

模仿示例

我当初阐明如何应用 JAGS 来模仿来自具备恒定生存和从新捕捉概率的模型的数据。我假如读者相熟这个模型及其作为状态空间模型的公式。

让咱们模仿一下!

# 恒定的生存和从新捕捉概率for (i in 1:nd){   for (t in f:(on-1)){#概率for (i in 1:nid){   # 定义埋伏状态和第一次捕捉时的察看值   z\[i,f\[i\] <- 1   mu2\[i,1\] <- 1 * z\[i,f\[i\] # 在第一次捕捉时检测为1("以第一次捕捉为条件")。   # 而后解决当前的状况   for (t in (f\[i\]+1):non){      # 状态过程      mu1\[i,t\] <- phi * z       # 察看过程      mu2\[i,t\] <- p * z

让咱们为参数抉择一些值并将它们存储在数据列表中:

# 用于模仿的参数 n = 100 # 个体的数量meanhi <- 0.8 # 存活率meap <- 0.6 # 重捕率data<-list

当初运行 JAGS

out

格式化输入:

as.mcmc(out)

head(dat)

我只监测了检测和非检测,但也能够取得状态的模仿值,即集体在每种状况下是生是死。你只须要批改对JAGS 的调用 monitor=c("y","x") 并相应地批改输入。

当初咱们将 Cormack-Jolly-Seber (CJS) 模型拟合到咱们刚刚模仿的数据中,假如参数不变:

# 倾向性和束缚for (i in 1:nd){   for (t in f\[i\]:(nn-1)){mehi ~ dunif(0, 1) # 均匀生存率的先验值Me ~ dunif(0, 1) # 均匀重捕的先验值# 概率for (i in 1:nd){   # 定义第一次捕捉时的埋伏状态   z\[i\]\] <- 1   for (t in (f\[i\]+1):nions){      # 状态过程      z\[i,t\] ~ dbern(mu1\[i,t\])      # 察看过程      y\[i,t\] ~ dbern(mu2\[i,t\])

筹备数据:

# 标记的场合的向量gerst <- function(x) min(which(x!=0))# 数据jagta
# 初始值for (i in 1:dim\]){       min(which(ch\[i,\]==1))      max(which(ch\[i,\]==1))function(){list(meaphi, mep , z ) }

咱们想对生存和从新捕捉的概率进行推断:

规范 MCMC 设置:

ni <- 10000

筹备运行 JAGS

# 从R中调用JAGSjags(nin = nb, woy = getwd() )

总结后验并与咱们用来模仿数据的值进行比拟:

print(cj3)

十分靠近!

跟踪图

trplot

后验分布图

denplot


最受欢迎的见解

1.matlab应用贝叶斯优化的深度学习

2.matlab贝叶斯隐马尔可夫hmm模型实现

3.R语言Gibbs抽样的贝叶斯简略线性回归仿真

4.R语言中的block Gibbs吉布斯采样贝叶斯多元线性回归

5.R语言中的Stan概率编程MCMC采样的贝叶斯模型

6.Python用PyMC3实现贝叶斯线性回归模型

7.R语言应用贝叶斯 层次模型进行空间数据分析

8.R语言随机搜寻变量抉择SSVS预计贝叶斯向量自回归(BVAR)模型

9.matlab贝叶斯隐马尔可夫hmm模型实现