RNG impact on MCMC [or lack thereof]
[This article was first published on R – Xi'an's Og, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.
Following the talk at MCM 2017 about the strange impact of the random generator on the outcome of an MCMC generator, I tried in Montréal airport the following code on the banana target of Haario et al. (1999), copied from Soetaert and Laine and using the MCMC function of the FME package:
library(FME) Banana <- function (x1, x2) { return(x2 - (x1^2+1)) } pmultinorm <- function(vec, mean, Cov) { diff <- vec - mean ex <- -0.5*t(diff) %*% solve(Cov) %*% diff rdet <- sqrt(det(Cov)) power <- -length(diff)*0.5 return((2.*pi)^power / rdet * exp(ex)) } BananaSS <- function (p) { P <- c(p[1], Banana(p[1], p[2])) Cov <- matrix(nr = 2, data = c(1, 0.9, 0.9, 1)) N=1e3 ejd=matrix(0,4,N) RNGkind("Mars") for (t in 1:N){ MCMC <- modMCMC(f = BananaSS, p = c(0, 0.7), jump = diag(nrow = 2, x = 5), niter = 1e3) ejd[1,t]=mean((MCMC$pars[-1,2]-MCMC$pars[1,2])^2)}
since this divergence from the initial condition seemed to reflect the experiment of the speaker at MCM 2017. Unsurprisingly, no difference came from using the different RNGs in R (which may fail to contain those incriminated by the study)…
Filed under: Books, R, Statistics, Travel, University life Tagged: Donald Knuth, George Marsaglia, GNU C library, MCM 2017, Montréal, R, random number generator, Super-duper
![](https://i0.wp.com/pixel.wp.com/b.gif?resize=1%2C1&ssl=1)
To leave a comment for the author, please follow the link and comment on their blog: R – Xi'an's Og.
R-bloggers.com offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.