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

From: Paul Gilbert <pgilbert_at_bank-banque-canada.ca>
Date: Thu 09 Jun 2005 - 01:19:38 EST

Beware that your easy trick will give you the same result every time you run it. You need a better scheme if you actually intend to get a new experiment each time you run it.

Paul

Gorjanc Gregor wrote:

> 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


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:38:06 2005

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