Replicate results with inla.posterior.sample

[Francisco Fabuel]

Hi all.

I would like to replicate the same results with inla.posterior.sample, but a code like this doesn't seem to work:

set.seed(1)

datos <- data.frame(x = 1:10)

datos$y <- 2 * datos$x + rnorm(10)

m <- inla(y ~ x, data = datos, control.compute=list(config = TRUE))

set.seed(1234)

l1 <- inla.posterior.sample(5, m)

set.seed(1234)

l2 <- inla.posterior.sample(5, m)

s1 <- sapply(l1, "[[", "latent")

s2 <- sapply(l2, "[[", "latent")

all.equal(s1, s2)

Thank you.

Francisco

[Håvard Rue]

The code also use an external C-program, for which random generator is
not aligned with the internal one in R. It would be possible to
regenerate the results using two rng seed values, and I can add that
feature if requested. 


H

On Sun, 2016-04-24 at 00:10 -0700, Francisco Fabuel wrote:
>
> set.seed(1)
> datos <- data.frame(x = 1:10)
> datos$y <- 2 * datos$x + rnorm(10)
>
> m <- inla(y ~ x, data = datos, control.compute=list(config = TRUE))
>
> set.seed(1234)
> l1 <- inla.posterior.sample(5, m)
>
> set.seed(1234)
> l2 <- inla.posterior.sample(5, m)
>
> s1 <- sapply(l1, "[[", "latent")
> s2 <- sapply(l2, "[[", "latent")

--
Håvard Rue
Department of Mathematical Sciences
Norwegian University of Science and Technology
N-7491 Trondheim, Norway
Voice: +47-7359-3533 URL : http://www.math.ntnu.no/~hrue
Mobile: +47-9260-0021 Email: havar...@math.ntnu.no

R-INLA: www.r-inla.org

[Francisco Fabuel]

I think there are situations where that possibility would be essential, for example, official statistics and replication of studies. So maybe you could consider to include it in future versions.

Thank you in advance.
Francisco

[Håvard Rue]

On Sun, 2016-04-24 at 02:58 -0700, Francisco Fabuel wrote:
> I think there are situations where that possibility would be
> essential, for example, official statistics and replication of
> studies. So maybe you could consider to include it in future
> versions.

I added this feature. It will appear in the new testing-version in a
few minutes. 

Thanks!
H

**********************************************

changeset:   3236:5b5d0435523a
tag:         tip
user:        Havard Rue <havar...@math.ntnu.no>
date:        Sun Apr 24 13:29:31 2016 +0200
files:       rinla/R/posterior.sample.R rinla/man/posterior.sample.Rd
description:
Added argument seed to inla.posterior.sample, so results can be
reproduced. Thanks to FF.

[Francisco Fabuel]

I've run the following code (R version 3.2.2, Platform: x86_64-w64-mingw32/x64 (64-bit)) and the samples are different.

set.seed(1)

datos <- data.frame(x = 1:10)

datos$y <- 2 * datos$x + rnorm(10)

m <- inla(y ~ x, data = datos, control.compute=list(config = TRUE))

germen <- 1234

l1 <- inla.posterior.sample(5, m, seed = germen)

l2 <- inla.posterior.sample(5, m, seed = germen)

s1 <- sapply(l1, "[[", "latent")

s2 <- sapply(l2, "[[", "latent")

all.equal(s1, s2)

[1] "Mean relative difference: 0.006738508"

First three columns of s1 and s2 are the same but something wrong happens with the other two.

Francisco

[Håvard Rue]

On Sun, 2016-04-24 at 23:09 -0700, Francisco Fabuel wrote:
> I've run the following code (R version 3.2.2, Platform: x86_64-w64-
> mingw32/x64 (64-bit)) and the samples are different.
>
> set.seed(1)
> datos <- data.frame(x = 1:10)
> datos$y <- 2 * datos$x + rnorm(10)
>
> m <- inla(y ~ x, data = datos, control.compute=list(config = TRUE))
>
> germen <- 1234
> l1 <- inla.posterior.sample(5, m, seed = germen)
> l2 <- inla.posterior.sample(5, m, seed = germen)
>
> s1 <- sapply(l1, "[[", "latent")
> s2 <- sapply(l2, "[[", "latent")
>
> all.equal(s1, s2)
> [1] "Mean relative difference: 0.006738508"
>
> First three columns of s1 and s2 are the same but something wrong
> happens with the other two.

see the example in ?inla.posterior.sample.

you need to control 2 seeds, in R and in `INLA' (which is a C-program).

do

> set.seed(1)
> datos <- data.frame(x = 1:10)
> datos$y <- 2 * datos$x + rnorm(10)
> 
> m <- inla(y ~ x, data = datos, control.compute=list(config = TRUE))
> 

> set.seed(8123)

> germen <- 1234
> l1 <- inla.posterior.sample(5, m, seed = germen)

> set.seed(8123)

> l2 <- inla.posterior.sample(5, m, seed = germen)
> 
> s1 <- sapply(l1, "[[", "latent")
> s2 <- sapply(l2, "[[", "latent")
> 
> all.equal(s1, s2)

[1] TRUE