[R] Random seed problem in MCMC coupling of chains

[Paul Gilbert]

The tools in setRNG are intended for this kind of problem and I do use 
them regularly in much more complicated situations.  They help save all 
the information, in addition to the seed, that you need for reproducible 
simulations. Try

niter <- 3
  nchain <- 2
  for (i in 1:niter) { # iterations
    tmpSeed <- setRNG()
    for (j in 1:nchain) { # chains
      setRNG(tmpSeed)
      a <- runif(1)
      cat("iter:", i, "chain:", j, "runif:", a, "\n")
    }
  }


iter: 1 chain: 1 runif: 0.8160078
iter: 1 chain: 2 runif: 0.8160078
iter: 2 chain: 1 runif: 0.4909793
iter: 2 chain: 2 runif: 0.4909793
iter: 3 chain: 1 runif: 0.4425924
iter: 3 chain: 2 runif: 0.4425924

HTH,
Paul Gilbert

Gorjanc Gregor wrote:

> Hello!
> 
> I am performing coupling of chains in MCMC and I need the same value
> of seed for two chains. I will show demo of what I want:
> 
> R code, which might show my example is:
> niter <- 3
> nchain <- 2
> tmpSeed <- 123
> for (i in 1:niter) { # iterations
>   for (j in 1:nchain) { # chains
>     set.seed(tmpSeed)
>     a <- runif(1)
>     cat("iter:", i, "chain:", j, "runif:", a, "\n")
>     tmpSeed <- .Random.seed
>   }
> }
> 
> I get this:
> 
> iter: 1 chain: 1 runif: 0.43588
> iter: 1 chain: 2 runif: 0.43588
> iter: 2 chain: 1 runif: 0.43588
> iter: 2 chain: 2 runif: 0.43588
> iter: 3 chain: 1 runif: 0.43588
> iter: 3 chain: 2 runif: 0.43588
> 
> but I would like to get:
> 
> iter: 1 chain: 1 runif: 0.43588
> iter: 1 chain: 2 runif: 0.43588
> iter: 2 chain: 1 runif: 0.67676
> iter: 2 chain: 2 runif: 0.67676
> iter: 3 chain: 1 runif: 0.12368
> iter: 3 chain: 2 runif: 0.12368
> 
> Note that seed value is of course changing, but it is parallel
> between chains.
> 
> I am able to do only this, since setting seed at the beginning 
> of chain i.e iteration is not a problem, but I want an upper 
> scheme, since I compare chains and stop one if some condition is
> satisfied.
> 
> tmpSeed <- 123
> for (i in 1:nchain) { # chains
>   set.seed(tmpSeed)
>   for (j in 1:niter) { # iterations
>     a <- runif(1)
>     cat("iter:", j, "chain:", i, "runif:", a, "\n")
>   }
> }
> iter: 1 chain: 1 runif: 0.28758
> iter: 2 chain: 1 runif: 0.7883
> iter: 3 chain: 1 runif: 0.40898
> iter: 1 chain: 2 runif: 0.28758
> iter: 2 chain: 2 runif: 0.7883
> iter: 3 chain: 2 runif: 0.40898
> iter: 1 chain: 3 runif: 0.28758
> iter: 2 chain: 3 runif: 0.7883
> iter: 3 chain: 3 runif: 0.40898
> 
> I was looking in 'rlecuyer', 'rsprng' and 'setRNG', but did not find 
> anything usable for me.  From reading on http://sprng.cs.fsu.edu/ 
> 'rsprng' provides just opposite of what I want, 'rlecuyer' is a bit
> to technical for me, but I think it also doesn't give identical
> seed for parallels. 'setRNG', especially it's function 'getRNG()'
> looks nice but its arguments should have seed stored. How can one
> do that?
> 
> 
> Thanks in advance!
> 
> Lep pozdrav / With regards,
>     Gregor Gorjanc
> 
> ----------------------------------------------------------------------
> University of Ljubljana
> Biotechnical Faculty        URI: http://www.bfro.uni-lj.si/MR/ggorjan
> Zootechnical Department     mail: gregor.gorjanc <at> bfro.uni-lj.si
> Groblje 3                   tel: +386 (0)1 72 17 861
> SI-1230 Domzale             fax: +386 (0)1 72 17 888
> Slovenia, Europe
> ----------------------------------------------------------------------
> "One must learn by doing the thing; for though you think you know it,
>  you have no certainty until you try." Sophocles ~ 450 B.C.
> 
> ______________________________________________
> R-help at stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html