Setting random seeds in rjags on a cluster

Help
2014-02-28
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
    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.