0

我在 R 中运行一个 for 循环（作为我使用 R2jags 运行的模型的功率分析的一部分）。在某些时候，我想知道我的 MCMC 链是否已经收敛，如果没有，我想跳过循环的迭代。但是，我不想跳到下一次迭代，我希望循环再次以相同的迭代开始。我目前正在使用命令“下一步”，但这正在跳过迭代。如何告诉我的 for 循环进行额外的迭代？下面是整个代码，但基本上我担心的就是这一点：

  if(up1 > 1.1 | up2 > 1.1 | up3 > 1.1 | up4 > 1.1)
   next

这是整个代码：

G.vec <- rep(NA, 1000)
#-----model------

model1 <- function(){  

  # likelihood
  for (c in 1:16){
    D.diff[c] ~ dnorm(mu, tau)
  }

  # priors
  mu ~ dnorm(0, 10)
  tau <- 1/(sd*sd)
  sd ~ dunif(0, 10)
}

#---- for loop--------

for (i in 1:1000){


#--- data simulation ---------
# paired data
sim_bf <- rtnorm(16, mean = 0.4949, sd = 0.12, lower = 0.1)
sim_af <- rtnorm(16, mean = 0.3959, sd = 0.12, lower = 0.1)
diff = sim_bf - sim_af


#------setting up jags--------------

data <- list(D.diff = diff)
params <- c("mu", "tau", "sd")
inits <- function(){
  list(mu = rnorm(1),
       sd = rlnorm(1)) 
}

#-------run jags----------------

output <- jags(data=data, 
               # inits=inits, 
               parameters.to.save=params,
               n.iter=1000, 
               n.burn=100, 
               n.chains=2, 
               n.thin=1,
               model.file=model1,
               progress.bar = "gui")

#-------convergence checker----------------
output.mcmc <- as.mcmc(output)
x <- gelman.diag(output.mcmc)
up1 <- x$psrf[1,2]        # Approximate convergence is diagnosed when the upper limit is close to 1
up2 <- x$psrf[2,2] 
up3 <- x$psrf[3,2] 
up4 <- x$psrf[4,2] 


  if(up1 > 1.1 | up2 > 1.1 | up3 > 1.1 | up4 > 1.1)
   next                 

# one sided t-test
  lo = output$BUGSoutput$summary[2,4]
  G.vec[i] <- ifelse((lo < 0), 0, 1)
}

a <- table(G.vec)
G <- a[2]/1000
G
4

2 回答 2

0

使用 break 语句作为收敛标准。

将下面的行放在 for 循环中是没有意义的，因为它预计会在收敛时运行。

lo = output$BUGSoutput$summary[2,4] G.vec[i] <- ifelse((lo < 0), 0, 1)

将它们放在 for 循环之外。

于 2015-04-29T05:21:55.937 回答
0

代替

for(i in 1:1000){...
    ...
    if(...
        next

做

i <- 0
while(i <=1000){...
    i <- i+1
    ...
    if(...
        i <- i-1
        break
于 2015-04-29T05:30:44.647 回答