[R] Random seed problem in MCMC coupling of chains

[Duncan Murdoch]

On 6/8/2005 9:27 AM, 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.

set.seed takes an integer, but .Random.seed is a complicated vector. 
You need to play with .Random.seed directly, and move your setting of 
tmpSeed out of the inner loop, i.e.

 > niter <- 3
 > nchain <- 2
 > set.seed(123)
 > tmpSeed <- .Random.seed
 > for (i in 1:niter) { # iterations
+   for (j in 1:nchain) { # chains
+     .Random.seed <- tmpSeed
+     a <- runif(1)
+     cat("iter:", i, "chain:", j, "runif:", a, "\n")
+   }
+   tmpSeed <- .Random.seed
+ }
iter: 1 chain: 1 runif: 0.2875775
iter: 1 chain: 2 runif: 0.2875775
iter: 2 chain: 1 runif: 0.7883051
iter: 2 chain: 2 runif: 0.7883051
iter: 3 chain: 1 runif: 0.4089769
iter: 3 chain: 2 runif: 0.4089769

However, heed the warnings in ?set.seed:  with some generators 
.Random.seed *does not* contain the full state of the RNG.  As far as I 
know there is no way to obtain the full state.

Duncan Murdoch