RE: [R] Random seed problem in MCMC coupling of chains

From: Gorjanc Gregor <Gregor.Gorjanc_at_bfro.uni-lj.si>
Date: Thu 09 Jun 2005 - 01:09:16 EST


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@stats.uwo.ca] Sent: sre 2005-06-08 15:53
To: Gorjanc Gregor
Cc: r-help@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



R-help@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 Received on Thu Jun 09 01:15:26 2005

This archive was generated by hypermail 2.1.8 : Fri 03 Mar 2006 - 03:32:27 EST