[R] Random seed problem in MCMC coupling of chains

[Gorjanc Gregor]

Thanks to Duncan, Dimitris as well as James for answers. I'll provide
here also example from James, which seems to be the easiest of them 
all and was not posted to the list:

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

Note that seed is set with iteration counter. This is really neat and
simple. I am just wondering if this is OK from "RNG point of view". Can
someone comment on that?

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.
----------------------------------------------------------------------

-----Original Message-----
From: Duncan Murdoch [mailto:murdoch at stats.uwo.ca]
Sent: sre 2005-06-08 15:53
To: Gorjanc Gregor
Cc: r-help at stat.math.ethz.ch
Subject: Re: [R] Random seed problem in MCMC coupling of chains
 
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