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:[EMAIL PROTECTED]
Sent: sre 2005-06-08 15:53
To: Gorjanc Gregor
Cc: [email protected]
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
______________________________________________
[email protected] mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
______________________________________________
[email protected] mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
