RUNJAGS-набор семян без предварительного моделирования - PullRequest
Новая
       92

RUNJAGS-набор семян без предварительного моделирования

0 голосов
/ 10 сентября 2018

Я выбираю некоторые данные из нормального дистрибутива, используя runjags. У меня нет никаких предварительных настроек для параметров, которые я использовал для своих симуляций. Кажется, что runjages не использует аргумент для исправления семени: list(".RNG.name"="base::Super-Duper", ".RNG.seed"=1). Я изменил аргумент на list(muOfClustsim=rep(1, npop), ".RNG.name"="base::Super-Duper", ".RNG.seed"=1), но он тоже не работает. Есть ли способ исправить семена для таких моделей в джунглях?

Вот минимальный воспроизводимый пример: 

library(runjags)

npop=3
nrep=10
sdpop=7
sigma=5
seed=4

set.seed(seed)
N = npop*nrep # nb of observations

## Population identity of each individual used to sample genotypes but not used for common garden test
pop <- rep(1:npop, each=nrep)

muOfClustsim <- rnorm(npop, 0, sdpop) # vector of population means
(tausim <- 1/(sigma*sigma)) # precision of random individual error

# parameters are treated as data for the simulation step
data <- list(N=N, pop=pop, muOfClustsim=muOfClustsim, tausim=tausim)

## JAG model
txtstring <- "
data{
  # Likelihood:
  for (i in 1:N){
    ysim[i] ~ dnorm(eta[i], tausim) # tau is precision (1 / variance)
    eta[i] <- muOfClustsim[pop[i]]
  }

}

model{
fake <- 0
}
"
## Initial values with seed for reproducibility
initssim <- list(".RNG.name"="base::Super-Duper", ".RNG.seed"=1)
##initssim <- list(muOfClustsim=rep(1, npop), ".RNG.name"="base::Super-Duper", ".RNG.seed"=1)
## Simulate with jags
set.seed(seed)
out <- run.jags(txtstring, data = data, monitor=c("ysim"), sample=1, n.chains=1, inits=initssim, summarise=FALSE)
## reformat the outputs
(ysim1 <- coda::as.mcmc(out)[1:N])

set.seed(seed)
out <- run.jags(txtstring, data = data, monitor=c("ysim"), sample=1, n.chains=1, inits=initssim, summarise=FALSE)
## reformat the outputs
(ysim2 <- coda::as.mcmc(out)[1:N])

identical(ysim1, ysim2)

1 Ответ

0 голосов
/ 11 сентября 2018

Исходные значения .RNG.name / .RNG.seed / .RNG.state применяются к модели (или, более конкретно, к цепочке внутри модели) и не используются в блоке данных. Это означает, что (насколько я знаю) нет способа сделать какие-либо стохастические представления в блоке данных, воспроизводимых в JAGS <= 4.3. Это то, что можно добавить в будущую версию JAGS, но я боюсь, что это будет низкий приоритет, так как всегда возможно (и обычно лучше) моделировать данные в R перед передачей их в JAGS. </p>

В вашем случае ответом (при условии, что вы хотите использовать JAGS) является симуляция в блоке модели, а не в блоке данных: 

txtstring <- "
model{
  # Likelihood:
  for (i in 1:N){
    ysim[i] ~ dnorm(eta[i], tausim) # tau is precision (1 / variance)
    eta[i] <- muOfClustsim[pop[i]]
  }
}
"

Остальная часть вашего кода запускается как ожидалось. Стоит отметить, что моделирование данных - это задача, обычно лучше подходящая для R, чем JAGS, но я предполагаю, что есть конкретная причина, по которой вы хотите использовать JAGS в этом случае ...

Мэтт

§ Хотя обычно не следует ожидать строгого равенства чисел, например: 

    identical(0.1+0.2, 0.3)

Но: 

    abs(0.3 - (0.1+0.2)) < sqrt(.Machine$double.eps)

Или еще лучше: 

    isTRUE(all.equal(0.1+0.2, 0.3))

Это стоит посмотреть: https://www.youtube.com/watch?v=3Bu7QUxzIbA&t=1s

...