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.
I am able to do only this, since setting seed at the beginning
of chain i.e iteration is not a problem, but I want an upper
scheme, since I compare chains and stop one if some condition is
satisfied.
tmpSeed <- 123
for (i in 1:nchain) { # chains
set.seed(tmpSeed)
for (j in 1:niter) { # iterations
a <- runif(1)
cat("iter:", j, "chain:", i, "runif:", a,
"\n")
}
}
iter: 1 chain: 1 runif: 0.28758
iter: 2 chain: 1 runif: 0.7883
iter: 3 chain: 1 runif: 0.40898
iter: 1 chain: 2 runif: 0.28758
iter: 2 chain: 2 runif: 0.7883
iter: 3 chain: 2 runif: 0.40898
iter: 1 chain: 3 runif: 0.28758
iter: 2 chain: 3 runif: 0.7883
iter: 3 chain: 3 runif: 0.40898
I was looking in 'rlecuyer', 'rsprng' and 'setRNG', but
did not find
anything usable for me. From reading on http://sprng.cs.fsu.edu/
'rsprng' provides just opposite of what I want, 'rlecuyer' is a
bit
to technical for me, but I think it also doesn't give identical
seed for parallels. 'setRNG', especially it's function
'getRNG()'
looks nice but its arguments should have seed stored. How can one
do that?
Thanks in advance!
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.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
Dimitris Rizopoulos
2005-Jun-08 13:53 UTC
[R] Random seed problem in MCMC coupling of chains
do you want something like this:
niter <- 3
nchain <- 2
rs <- sample(500, niter, TRUE)
for (i in 1:niter) { # iterations
for (j in 1:nchain) { # chains
set.seed(rs[i])
a <- runif(1)
cat("iter:", i, "chain:", j, "runif:", a,
"\n")
}
}
I hope it helps.
Best,
Dimitris
----
Dimitris Rizopoulos
Ph.D. Student
Biostatistical Centre
School of Public Health
Catholic University of Leuven
Address: Kapucijnenvoer 35, Leuven, Belgium
Tel: +32/16/336899
Fax: +32/16/337015
Web: http://www.med.kuleuven.ac.be/biostat/
http://www.student.kuleuven.ac.be/~m0390867/dimitris.htm
----- Original Message -----
From: "Gorjanc Gregor" <Gregor.Gorjanc at bfro.uni-lj.si>
To: <r-help at stat.math.ethz.ch>
Sent: Wednesday, June 08, 2005 3:27 PM
Subject: [R] Random seed problem in MCMC coupling of chains
> 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.
>
> I am able to do only this, since setting seed at the beginning
> of chain i.e iteration is not a problem, but I want an upper
> scheme, since I compare chains and stop one if some condition is
> satisfied.
>
> tmpSeed <- 123
> for (i in 1:nchain) { # chains
> set.seed(tmpSeed)
> for (j in 1:niter) { # iterations
> a <- runif(1)
> cat("iter:", j, "chain:", i, "runif:", a,
"\n")
> }
> }
> iter: 1 chain: 1 runif: 0.28758
> iter: 2 chain: 1 runif: 0.7883
> iter: 3 chain: 1 runif: 0.40898
> iter: 1 chain: 2 runif: 0.28758
> iter: 2 chain: 2 runif: 0.7883
> iter: 3 chain: 2 runif: 0.40898
> iter: 1 chain: 3 runif: 0.28758
> iter: 2 chain: 3 runif: 0.7883
> iter: 3 chain: 3 runif: 0.40898
>
> I was looking in 'rlecuyer', 'rsprng' and 'setRNG',
but did not find
> anything usable for me. From reading on http://sprng.cs.fsu.edu/
> 'rsprng' provides just opposite of what I want, 'rlecuyer'
is a bit
> to technical for me, but I think it also doesn't give identical
> seed for parallels. 'setRNG', especially it's function
'getRNG()'
> looks nice but its arguments should have seed stored. How can one
> do that?
>
>
> Thanks in advance!
>
> 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.
>
> ______________________________________________
> R-help at 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
>
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
The tools in setRNG are intended for this kind of problem and I do use
them regularly in much more complicated situations. They help save all
the information, in addition to the seed, that you need for reproducible
simulations. Try
niter <- 3
nchain <- 2
for (i in 1:niter) { # iterations
tmpSeed <- setRNG()
for (j in 1:nchain) { # chains
setRNG(tmpSeed)
a <- runif(1)
cat("iter:", i, "chain:", j, "runif:", a,
"\n")
}
}
iter: 1 chain: 1 runif: 0.8160078
iter: 1 chain: 2 runif: 0.8160078
iter: 2 chain: 1 runif: 0.4909793
iter: 2 chain: 2 runif: 0.4909793
iter: 3 chain: 1 runif: 0.4425924
iter: 3 chain: 2 runif: 0.4425924
HTH,
Paul Gilbert
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.
>
> I am able to do only this, since setting seed at the beginning
> of chain i.e iteration is not a problem, but I want an upper
> scheme, since I compare chains and stop one if some condition is
> satisfied.
>
> tmpSeed <- 123
> for (i in 1:nchain) { # chains
> set.seed(tmpSeed)
> for (j in 1:niter) { # iterations
> a <- runif(1)
> cat("iter:", j, "chain:", i, "runif:", a,
"\n")
> }
> }
> iter: 1 chain: 1 runif: 0.28758
> iter: 2 chain: 1 runif: 0.7883
> iter: 3 chain: 1 runif: 0.40898
> iter: 1 chain: 2 runif: 0.28758
> iter: 2 chain: 2 runif: 0.7883
> iter: 3 chain: 2 runif: 0.40898
> iter: 1 chain: 3 runif: 0.28758
> iter: 2 chain: 3 runif: 0.7883
> iter: 3 chain: 3 runif: 0.40898
>
> I was looking in 'rlecuyer', 'rsprng' and 'setRNG',
but did not find
> anything usable for me. From reading on http://sprng.cs.fsu.edu/
> 'rsprng' provides just opposite of what I want, 'rlecuyer'
is a bit
> to technical for me, but I think it also doesn't give identical
> seed for parallels. 'setRNG', especially it's function
'getRNG()'
> looks nice but its arguments should have seed stored. How can one
> do that?
>
>
> Thanks in advance!
>
> 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.
>
> ______________________________________________
> R-help at 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
Thanks to Paul and Gabor for additional tips/examples. Actually, I find
Pauls suggestion with setRNG also nice and is exactly what I wanted.
Paul, if I understand this correctly, your suggestion with setRNG does not
alter "RNG flow", it just takes care that chains really have equal
seeds.
I remember that I have read somewhere that destroying "RNG flow over and
over to get real randomness" is not a good idea. Can someone confirm this?
niter <- 3
nchain <- 2
for (i in 1:niter) { # iterations
tmpSeed <- setRNG()
for (j in 1:nchain) { # chains
setRNG(tmpSeed)
a <- runif(1)
cat("iter:", i, "chain:", j, "runif:", a,
"\n")
}
}
iter: 1 chain: 1 runif: 0.8160078
iter: 1 chain: 2 runif: 0.8160078
iter: 2 chain: 1 runif: 0.4909793
iter: 2 chain: 2 runif: 0.4909793
iter: 3 chain: 1 runif: 0.4425924
iter: 3 chain: 2 runif: 0.4425924
[... removed other stuff ...]
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.
Gabor Grothendieck
2005-Jun-08 18:42 UTC
[R] Random seed problem in MCMC coupling of chains
Here is a small variation. We define a list to hold
the last seed for each chain. Each time we enter the simulation
for a chain we use that seed and each time we exit we update it.
The loop becomes simpler since the setup is all done prior
to looping and everything else is done in the inner loop.
Note that a double loop with nothing between the first and second
for is really like a single loop over the i,j pairs so its presumably
easier to understand.
library(setRNG)
set.seed(123)
niter <- 3; nchain <- 2
chain <- lapply(1:nchain, function(x) setRNG())
for(i in 1:niter)
for(j in 1:nchain) {
setRNG(chain[[j]]) # get seed
a <- runif(1)
cat("iter:", i, "chain:", j, "runif:", a,
"\n")
chain[[j]] <- setRNG() # save seed
}
On 6/8/05, Gorjanc Gregor <Gregor.Gorjanc at bfro.uni-lj.si>
wrote:> Thanks to Paul and Gabor for additional tips/examples. Actually, I find
> Pauls suggestion with setRNG also nice and is exactly what I wanted.
> Paul, if I understand this correctly, your suggestion with setRNG does not
> alter "RNG flow", it just takes care that chains really have
equal seeds.
> I remember that I have read somewhere that destroying "RNG flow over
and
> over to get real randomness" is not a good idea. Can someone confirm
this?
>
> niter <- 3
> nchain <- 2
> for (i in 1:niter) { # iterations
> tmpSeed <- setRNG()
> for (j in 1:nchain) { # chains
> setRNG(tmpSeed)
> a <- runif(1)
> cat("iter:", i, "chain:", j, "runif:", a,
"\n")
> }
> }
>
> iter: 1 chain: 1 runif: 0.8160078
> iter: 1 chain: 2 runif: 0.8160078
> iter: 2 chain: 1 runif: 0.4909793
> iter: 2 chain: 2 runif: 0.4909793
> iter: 3 chain: 1 runif: 0.4425924
> iter: 3 chain: 2 runif: 0.4425924
>
> [... removed other stuff ...]
>
> 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.
> ----------------------------------------------------------------------
>
>
>
>
>
>
>
>