# 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

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.

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.

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