gpt4 book ai didi

r - 对各种迭代应用贝叶斯模型 (JAGS)

转载 作者:行者123 更新时间:2023-12-03 16:01:58 26 4
gpt4 key购买 nike

考虑以下数据框:

set.seed(5678)
sub_df<- data.frame(clustersize= rep(1, 4),
lepsp= c("A", "B", "C", "D"),
dens= round(runif(4, c(0, 1)), 3),
db= sample(1:10, 4, replace=TRUE))
假设我想运行以下贝叶斯线性模型,它返回 samples , 一个 mc.array目的:
library("rjags")
library("coda")
dataForJags <- list(dens=sub_df$dens, db=sub_df$db, N=length(sub_df$dens))


model<-"model{
for(i in 1:N){
dens[i] ~ dnorm(mu[i], tau)
# identity
mu[i] <- int + beta1*db[i]
}
tau ~ dgamma(0.1,0.1)
int ~ dnorm(0, 0.001)
beta1 ~ dnorm(0, 0.001)
}"

##compile

mod1 <- jags.model(textConnection(model),data= dataForJags,n.chains=2)

##samples returns a list of mcarray objects

samples<-jags.samples(model= mod1,variable.names=c("beta1",
"int","mu","tau"),n.iter=100000)
鉴于 samples$beta1[,,]表示来自 jags 模型参数的后验分布的随机样本,然后总结一下,我的下一步是计算后验分布的均值和 95% 可信区间。所以我会使用:
coeff_output<- round(quantile(samples$beta1[,,],probs=c(0.5,0.025,0.975)),3)
现在,假设我的实际数据框有多个级别 clustersize .
set.seed(5672)
df<- data.frame(clustersize= c(rep(1, 4), rep(2,4), rep(3, 3)),
lepsp= c("A", "B", "C", "D", "B", "C", "D", "E", "A", "D", "F"),
dens= round(runif(11, c(0, 1)), 3),
db= sample(1:10, 11, replace=TRUE))
我将如何为 clustersize 的每个级别运行此模型?单独并使用 forloop 将输出编译为单个结果数据帧或 apply功能?对于 clustersize 的每个级别, 结果 mc.array对象 samples应该输出到 result_listcoeff_output应该输出到数据框 result_coeff .
下面我计算每个 clustersize 的输出分别生成预期结果列表和数据框。
 #clustersize==1
sub_df1<- data.frame(clustersize= rep(1, 4),
lepsp= c("A", "B", "C", "D"),
dens= round(runif(4, c(0, 1)), 3),
db= sample(1:10, 4, replace=TRUE))

dataForJags <- list(dens=sub_df$dens, db=sub_df$db, N=length(sub_df$dens))
model<-"model{
for(i in 1:N){
dens[i] ~ dnorm(mu[i], tau)
mu[i] <- int + beta1*db[i]
}
tau ~ dgamma(0.1,0.1)
int ~ dnorm(0, 0.001)
beta1 ~ dnorm(0, 0.001)
}"

mod1 <- jags.model(textConnection(model),data= dataForJags,n.chains=2)

samples1<-jags.samples(model= mod1,variable.names=c("beta1",
"int","mu","tau"),n.iter=100000)

coeff_output1<-
data.frame(as.list(round(quantile(samples1$beta1[,,],probs=c(0.5,0.025,0.975)),3)))

#clustersize==2
sub_df2<- data.frame(clustersize= rep(2,4),
lepsp= c( "B", "C", "D", "E"),
dens= round(runif(4, c(0, 1)), 3),
db= sample(1:10, 4, replace=TRUE))
dataForJags <- list(dens=sub_df$dens, db=sub_df$db, N=length(sub_df$dens))
model<-"model{
for(i in 1:N){
dens[i] ~ dnorm(mu[i], tau)
mu[i] <- int + beta1*db[i]
}
tau ~ dgamma(0.1,0.1)
int ~ dnorm(0, 0.001)
beta1 ~ dnorm(0, 0.001)
}"

mod1 <- jags.model(textConnection(model),data= dataForJags,n.chains=2)

samples2<-jags.samples(model= mod1,variable.names=c("beta1",
"int","mu","tau"),n.iter=100000)

coeff_output2<-
data.frame(as.list(round(quantile(samples2$beta1[,,],probs=c(0.5,0.025,0.975)),3)))

#clustersize==3
sub_df3<- data.frame(clustersize= rep(3, 3),
lepsp= c("A", "D", "F"),
dens= round(runif(3, c(0, 1)), 3),
db= sample(1:10, 3, replace=TRUE))
dataForJags <- list(dens=sub_df$dens, db=sub_df$db, N=length(sub_df$dens))
model<-"model{
for(i in 1:N){
dens[i] ~ dnorm(mu[i], tau)
mu[i] <- int + beta1*db[i]
}
tau ~ dgamma(0.1,0.1)
int ~ dnorm(0, 0.001)
beta1 ~ dnorm(0, 0.001)
}"

mod1 <- jags.model(textConnection(model),data= dataForJags,n.chains=2)

samples3<-jags.samples(model= mod1,variable.names=c("beta1",
"int","mu","tau"),n.iter=100000)

coeff_output3<-
data.frame(as.list(round(quantile(samples3$beta1[,,],probs=c(0.5,0.025,0.975)),3)))
期望的最终输出:
result_list<- list(samples1, samples2, samples3)

result_coeff<-rbind(coeff_output1, coeff_output2, coeff_output3)
这是实际数据框的链接。该解决方案应该能够处理集群大小高达 600 的大型数据帧。
download.file("https://drive.google.com/file/d/1ZYIQtb_QHbYsInDGkta-5P2EJrFRDf22/view?usp=sharing",temp)

最佳答案

这里有几个问题需要考虑,这些问题是由您尝试做的事情的规模引起的。您正在创建超过 550 个不同的 jags.sample每个对象都有 100000 次迭代,然后尝试将它们全部存储在一个列表中。在大多数机器上,这会导致内存问题:输出太大了。
我们至少有两种方法可以解决这个问题:

  • 采取措施尽可能减少我们输入数据的内存使用量。
  • 调整我们的 JAGS 输出,这样它就不会从每个链中保存这么多的迭代。

  • 我对您的代码进行了许多修改,以使其能够与您的实际数据集一起使用。
    创建输入数据:
    在您的原始代码中,clustersize 和 db 都具有数据类型 numeric ,即使它们只需要是整数。 numeric类型占用 8 个字节,而 integer type 只占用 4 个字节。如果我们将这两列强制为 integer类型,我们实际上可以将下一步数据帧列表的内存大小减少大约 30%。
    library("tidyverse")

    #### Load Raw Data ####
    df <- read_csv("example.csv") %>%
    select(-1) %>%
    mutate(clustersize = as.integer(clustersize),
    db = as.integer(db))
    初始 JAGS 调整
    您为每个链使用了太多的迭代; niter = 100000 非常高。您还应该使用 n.burn 指定老化期。 , 使用 n.adapt 的适应期,以及使用 thin 的细化参数.细化参数在这里特别重要——这直接减少了我们从每个链中节省的迭代次数。细化参数为 50 意味着我们只保存每 50 个结果。
    有一些事后方法可以选择细化参数、老化和适应期,但这种讨论超出了 SO 的范围。有关所有这些论点的作用的一些基本信息,这里有一个很好的答案: https://stackoverflow.com/a/38875637/9598813 .目前,我提供了允许此代码在整个数据集上运行的值,但我建议您仔细选择用于最终分析的值。
    使用 tidybayes
    以下解决方案使用 tidybayes包裹。这提供了一个干净的输出,并允许我们巧妙地将所有系数汇总行绑定(bind)到单个数据帧中。请注意,我们使用 coda.samples()而不是 jags.samples() ,因为这提供了一个更通用的 MCMC 对象,我们可以将其传递给 spread_draws() .我们还使用 dplyr::group_split()这比 split() 的计算效率略高.
    library("rjags")
    library("coda")
    library("tidybayes")

    set.seed(5672)
    result <- df %>% group_split(clustersize) %>% map(~{

    dataForJags <- list(dens=.x$dens, db=.x$db, N=length(.x$dens))

    # Declare model structure
    mod1 <- jags.model(textConnection(model),
    data=dataForJags,
    n.chains=2)

    # samples returns a list of mcmc objects
    samples<-coda.samples(model=mod1,
    variable.names=c("beta1","int","mu","tau"),
    n.burn=10000,
    n.adapt=5000,
    n.iter=25000,
    thin=50
    )
    # Extract individual draws
    samp <- spread_draws(samples, beta1)

    # Summarize 95% credible intervals
    coeff_output <- spread_draws(samples, beta1) %>%
    median_qi(beta1)

    list(samples = samp, coeff_output = coeff_output)
    }) %>% transpose()

    # List of sample objects
    result$samples
    # Dataframe of coefficient estimates and 95% credible intervals
    result_coeff <- bind_rows(result$coeff_output, .id = "clustersize")

    关于r - 对各种迭代应用贝叶斯模型 (JAGS),我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/62648243/

    26 4 0
    Copyright 2021 - 2024 cfsdn All Rights Reserved 蜀ICP备2022000587号
    广告合作:1813099741@qq.com 6ren.com