标签：Seber MCMC 模型 Jolly 参数 JAGS 贝叶斯 数据 模拟

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

原文出处：拓端数据部落公众号

本文，我通过两个种群生态学家可能感兴趣的例子来说明使用“JAGS”来模拟数据：首先是线性回归，其次是估计动物存活率（公式化为状态空间模型）。

最近，我一直在努力模拟来自复杂分层模型的数据。我现在正在使用 JAGS。

模拟数据 JAGS 很方便，因为你可以使用（几乎）相同的代码进行模拟和推理，并且你可以在相同的环境（即JAGS）中进行模拟研究（偏差、精度、区间覆盖 ）。

线性回归示例

我们首先加载本教程所需的包：

1.   library(R2jags)
2.    
3.    

然后直接切入正题，让我们从线性回归模型生成数据。使用一个 data 块，并将参数作为数据传递。

1.   data{
2.   # 似然函数:
3.   for (i in 1:N){
4.   y[i] ~ # tau是精度（1/方差）。
5.    
6.   }
7.    

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

我们为用作数据的模型参数选择一些值：

1.    
2.   # 模拟的参数
3.   N # 样本
4.   x <- 1:N # 预测因子
5.   alpha # 截距
6.   beta # 斜率
7.   sigma# 残差sd
8.   1/(sigma*sigma) # 精度
9.   # 在模拟步骤中，参数被当作数据处理
10.    

现在运行 JAGS; 请注意，我们监控因变量而不是参数，就像我们在进行标准推理时所做的那样：

1.   # 运行结果
2.   out

输出有点乱，需要适当格式化：

1.   # 重新格式化输出
2.   mcmc(out)

dim

dat

现在让我们将我们用来模拟的模型拟合到我们刚刚生成的数据中。不再赘述，假设读者熟悉 JAGS 线性回归。

1.    
2.   # 用BUGS语言指定模型
3.   model <-
4.    
5.   for (i in 1:N){
6.   y[i] ~ dnorm(mu[i], tau) # tau是精度(1/方差)
7.    
8.    
9.   alpha 截距
10.   beta # 斜率
11.   sigma # 标准差
12.    
13.    
14.   # 数据
15.   dta <- list(y = dt, N = length(at), x = x)
16.    
17.   # 初始值
18.   inits
19.    
20.    
21.   # MCMC设置
22.   ni <- 10000
23.    
24.    
25.   # 从R中调用JAGS
26.   jags()
27.    

让我们看看结果并与我们用来模拟数据的参数进行比较（见上文）：

1.   # 总结后验
2.   print(res)

检查收敛：

1.   # 追踪图
2.   plot(res)

绘制回归参数和残差标准差的后验分布：

1.   # 后验分布
2.   plot(res)

模拟示例

我现在说明如何使用 JAGS 来模拟来自具有恒定生存和重新捕获概率的模型的数据。我假设读者熟悉这个模型及其作为状态空间模型的公式。

让我们模拟一下！

1.    
2.   # 恒定的生存和重新捕获概率
3.   for (i in 1:nd){
4.   for (t in f:(on-1)){
5.    
6.   #概率
7.   for (i in 1:nid){
8.   # 定义潜伏状态和第一次捕获时的观察值
9.   z[i,f[i] <- 1
10.   mu2[i,1] <- 1 * z[i,f[i] # 在第一次捕获时检测为1（"以第一次捕获为条件"）。
11.   # 然后处理以后的情况
12.   for (t in (f[i]+1):non){
13.   # 状态进程
14.   mu1[i,t] <- phi * z
15.   # 观察过程
16.   mu2[i,t] <- p * z
17.    
18.    

让我们为参数选择一些值并将它们存储在数据列表中：

1.    
2.   # 用于模拟的参数
3.   n = 100 # 个体的数量
4.   meanhi <- 0.8 # 存活率
5.   meap <- 0.6 # 重捕率
6.   data<-list

现在运行 JAGS：

out 

格式化输出：

as.mcmc(out)

head(dat)

我只监测了检测和非检测，但也可以获得状态的模拟值，即个人在每种情况下是生是死。你只需要修改对JAGS 的调用 monitor=c("y","x") 并相应地修改输出。

现在我们将 Cormack-Jolly-Seber (CJS) 模型拟合到我们刚刚模拟的数据中，假设参数不变：

1.    
2.   # 倾向性和约束
3.   for (i in 1:nd){
4.   for (t in f[i]:(nn-1)){
5.    
6.    
7.   mehi ~ dunif(0, 1) # 平均生存率的先验值
8.   Me ~ dunif(0, 1) # 平均重捕的先验值
9.   # 概率
10.   for (i in 1:nd){
11.   # 定义第一次捕获时的潜伏状态
12.   z[i]] <- 1
13.   for (t in (f[i]+1):nions){
14.   # 状态过程
15.   z[i,t] ~ dbern(mu1[i,t])
16.   # 观察过程
17.   y[i,t] ~ dbern(mu2[i,t])

准备数据：

1.    
2.   # 标记的场合的向量
3.   gerst <- function(x) min(which(x!=0))
4.   # 数据
5.   jagta

1.    
2.   # 初始值
3.   for (i in 1:dim]){
4.   min(which(ch[i,]==1))
5.   max(which(ch[i,]==1))
6.    
7.   function(){list(meaphi, mep , z ) }
8.    

我们想对生存和重新捕获的概率进行推断：

标准 MCMC 设置：

ni <- 10000

准备运行 JAGS！

1.   # 从R中调用JAGS
2.   jags(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模型实现

标签：Seber,MCMC,模型,Jolly,参数,JAGS,贝叶斯,数据,模拟
来源： https://www.cnblogs.com/tecdat/p/15843753.html

本站声明： 1. iCode9 技术分享网（下文简称本站）提供的所有内容，仅供技术学习、探讨和分享；
2. 关于本站的所有留言、评论、转载及引用，纯属内容发起人的个人观点，与本站观点和立场无关；
3. 关于本站的所有言论和文字，纯属内容发起人的个人观点，与本站观点和立场无关；
4. 本站文章均是网友提供，不完全保证技术分享内容的完整性、准确性、时效性、风险性和版权归属；如您发现该文章侵犯了您的权益，可联系我们第一时间进行删除；
5. 本站为非盈利性的个人网站，所有内容不会用来进行牟利，也不会利用任何形式的广告来间接获益，纯粹是为了广大技术爱好者提供技术内容和技术思想的分享性交流网站。