R語(yǔ)言生態(tài)學(xué)JAGS模擬數(shù)據(jù)、線性回歸、Cormack-Jolly-Seber (CJS) 模型擬合MCMC 估計(jì)動(dòng)
原文鏈接：http://tecdat.cn/?p=24721??
原文出處：拓端數(shù)據(jù)部落公眾號(hào)
本文，我通過兩個(gè)種群生態(tài)學(xué)家可能感興趣的例子來(lái)說(shuō)明使用“JAGS”來(lái)模擬數(shù)據(jù)：首先是線性回歸，其次是估計(jì)動(dòng)物存活率（公式化為狀態(tài)空間模型）。
最近，我一直在努力模擬來(lái)自復(fù)雜分層模型的數(shù)據(jù)。我現(xiàn)在正在使用?JAGS。
模擬數(shù)據(jù)?JAGS?很方便，因?yàn)槟憧梢允褂茫◣缀酰┫嗤拇a進(jìn)行模擬和推理，并且你可以在相同的環(huán)境（即JAGS）中進(jìn)行模擬研究（偏差、精度、區(qū)間覆蓋?）。
線性回歸示例
我們首先加載本教程所需的包：
library(R2jags)
然后直接切入正題，讓我們從線性回歸模型生成數(shù)據(jù)。使用一個(gè)?data?塊，并將參數(shù)作為數(shù)據(jù)傳遞。
data{# 似然函數(shù):for (i in 1:N){y[i] ~ ?# tau是精度（1/方差）。}
這里，?alpha?和?beta?是截距和斜率、?tau?方差的精度或倒數(shù)、?y?因變量和?x?解釋變量。
我們?yōu)橛米鲾?shù)據(jù)的模型參數(shù)選擇一些值：
# 模擬的參數(shù)N ?# 樣本x <- 1:N # 預(yù)測(cè)因子alpha # 截距beta ?# 斜率sigma# 殘差sd1/(sigma*sigma) # 精度# 在模擬步驟中，參數(shù)被當(dāng)作數(shù)據(jù)處理
現(xiàn)在運(yùn)行?JAGS; 請(qǐng)注意，我們監(jiān)控因變量而不是參數(shù)，就像我們?cè)谶M(jìn)行標(biāo)準(zhǔn)推理時(shí)所做的那樣：
# 運(yùn)行結(jié)果out
輸出有點(diǎn)亂，需要適當(dāng)格式化：
# 重新格式化輸出mcmc(out)
dim
dat
現(xiàn)在讓我們將我們用來(lái)模擬的模型擬合到我們剛剛生成的數(shù)據(jù)中。不再贅述，假設(shè)讀者熟悉?JAGS?線性回歸。
# 用BUGS語(yǔ)言指定模型model <-for (i in 1:N){y[i] ~ dnorm(mu[i], tau) # tau是精度(1/方差)alpha ?截距beta # 斜率sigma ?# 標(biāo)準(zhǔn)差# 數(shù)據(jù)dta <- list(y = dt, N = length(at), x = x)# 初始值inits# MCMC設(shè)置ni <- 10000# 從R中調(diào)用JAGSjags()
讓我們看看結(jié)果并與我們用來(lái)模擬數(shù)據(jù)的參數(shù)進(jìn)行比較（見上文）：
# 總結(jié)后驗(yàn)print(res)
檢查收斂：
# 追蹤圖plot(res)
繪制回歸參數(shù)和殘差標(biāo)準(zhǔn)差的后驗(yàn)分布：
# 后驗(yàn)分布plot(res)
模擬示例
我現(xiàn)在說(shuō)明如何使用?JAGS?來(lái)模擬來(lái)自具有恒定生存和重新捕獲概率的模型的數(shù)據(jù)。我假設(shè)讀者熟悉這個(gè)模型及其作為狀態(tài)空間模型的公式。
讓我們模擬一下！
# 恒定的生存和重新捕獲概率for (i in 1:nd){for (t in f:(on-1)){#概率for (i in 1:nid){# 定義潛伏狀態(tài)和第一次捕獲時(shí)的觀察值z[i,f[i] <- 1mu2[i,1] <- 1 * z[i,f[i] # 在第一次捕獲時(shí)檢測(cè)為1（"以第一次捕獲為條件"）。# 然后處理以后的情況for (t in (f[i]+1):non){# 狀態(tài)進(jìn)程mu1[i,t] <- phi * z# 觀察過程mu2[i,t] <- p * z
讓我們?yōu)閰?shù)選擇一些值并將它們存儲(chǔ)在數(shù)據(jù)列表中：
# 用于模擬的參數(shù)n = 100 # 個(gè)體的數(shù)量meanhi <- 0.8 # 存活率meap <- 0.6 # 重捕率data<-list
現(xiàn)在運(yùn)行?JAGS：
out
格式化輸出：
as.mcmc(out)
head(dat)
我只監(jiān)測(cè)了檢測(cè)和非檢測(cè)，但也可以獲得狀態(tài)的模擬值，即個(gè)人在每種情況下是生是死。你只需要修改對(duì)JAGS?的調(diào)用?monitor=c("y","x")?并相應(yīng)地修改輸出。
現(xiàn)在我們將 Cormack-Jolly-Seber (CJS) 模型擬合到我們剛剛模擬的數(shù)據(jù)中，假設(shè)參數(shù)不變：
# 傾向性和約束for (i in 1:nd){for (t in f[i]:(nn-1)){mehi ~ dunif(0, 1) # 平均生存率的先驗(yàn)值Me ~ dunif(0, 1) # 平均重捕的先驗(yàn)值# 概率for (i in 1:nd){# 定義第一次捕獲時(shí)的潛伏狀態(tài)z[i]] <- 1for (t in (f[i]+1):nions){# 狀態(tài)過程z[i,t] ~ dbern(mu1[i,t])# 觀察過程y[i,t] ~ dbern(mu2[i,t])
準(zhǔn)備數(shù)據(jù)：
# 標(biāo)記的場(chǎng)合的向量gerst <- function(x) min(which(x!=0))# 數(shù)據(jù)jagta
# 初始值
for (i in 1:dim]){
min(which(ch[i,]==1))
max(which(ch[i,]==1))
function(){list(meaphi, mep , z ) }
我們想對(duì)生存和重新捕獲的概率進(jìn)行推斷：
標(biāo)準(zhǔn) MCMC 設(shè)置：
ni <- 10000
準(zhǔn)備運(yùn)行?JAGS！
# 從R中調(diào)用JAGSjags(nin = nb, woy = getwd() )
總結(jié)后驗(yàn)并與我們用來(lái)模擬數(shù)據(jù)的值進(jìn)行比較：
print(cj3)
非常接近！
跟蹤圖
trplot
后驗(yàn)分布圖
denplot
最受歡迎的見解
1.matlab使用貝葉斯優(yōu)化的深度學(xué)習(xí)
2.matlab貝葉斯隱馬爾可夫hmm模型實(shí)現(xiàn)
3.R語(yǔ)言Gibbs抽樣的貝葉斯簡(jiǎn)單線性回歸仿真
4.R語(yǔ)言中的block Gibbs吉布斯采樣貝葉斯多元線性回歸
5.R語(yǔ)言中的Stan概率編程MCMC采樣的貝葉斯模型
6.Python用PyMC3實(shí)現(xiàn)貝葉斯線性回歸模型
7.R語(yǔ)言使用貝葉斯 層次模型進(jìn)行空間數(shù)據(jù)分析
8.R語(yǔ)言隨機(jī)搜索變量選擇SSVS估計(jì)貝葉斯向量自回歸（BVAR）模型
9.matlab貝葉斯隱馬爾可夫hmm模型實(shí)現(xiàn)
