Torbjørn Ergon
2014-02-28
I running Monte-Carlo simulations involving JAGS run through rjags on a computer cluster. I understand that JAGS sets the random seeds from the current time. Hence I supply ".RNG.name" and ".RNG.seed" in the list of initial parameters, but this does not seem to work: For one of my models and data-subsets, where I expect that there is rather little information about the parameters in the data, I get lots of identical mcmc-chains. Moreover, the identical chains have always been started around the same time on different cpu's. Initial values and data are always different (but I think that some simulated data have very little information about the parameters in the model). I pass a unique seeds to R on each worker and confirmed that this works.
What I do is this:
Inits.1sp.1meth = function(z, psi, p, seed){ list( psi = invlogit(log((psi/(1-psi))*runif(1,0.8,1.2))), p = invlogit(log((p/(1-p))*runif(1,0.8,1.2))), z = z, .RNG.name = "base::Mersenne-Twister", .RNG.seed = seed ) } seeds = matrix(sample(1:1000000, 8*3), 8, 3) # I've provided different seeds for each worker and confirmed that these are indeed unique for each run. inits.spA.1meth = list( Inits.1sp.1meth(z[,1], psi[1], p2[1], seeds[7,1]), Inits.1sp.1meth(z[,1], psi[1], p2[1], seeds[7,2]), Inits.1sp.1meth(z[,1], psi[1], p2[1], seeds[7,3]) ) setup.spA.1meth = jags.model(file = paste(model.path, "1sp1meth.txt", sep = ""), data = DataSim.spA.1meth, inits = inits.spA.1meth, n.chains = 3, n.adapt=3000) jags.spA.1meth = coda.samples(setup.spA.1meth, variable.names = params.1meth, n.iter=3000, thin=1)
It does not seem to be possible to provide a seed for coda.samples(), so I guess this is passed on from the jags model object (or could the problem be that it is not?).
Any advice or hints would be deeply appreciated.
Best regards,
Torbjørn Ergon
Martyn Plummer
2014-02-28
To check that the initializing of the RNG is working correctly, use n.adapt=0
when you create the model
setup.spA.1meth = jags.model(file = paste(model.path, "1sp1meth.txt", sep = ""), data = DataSim.spA.1meth, inits = inits.spA.1meth, n.chains = 3, n.adapt=0)
Then immediately call
setup.spA.1meth$state(internal=TRUE)
You should get a reproducible state (although you will not be able to trace the translation from .RNG.seed
to .RNG.state
the latter is a deterministic function of the former).
The parallel.seeds()
function can be used to generate random seeds. See also the dclone
package, which has some facilities for running JAGS in parallel.