Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.
This special choice of a single proposal however simplifies the acceptance ratio (and keeps the overall proposal symmetric). The final acceptance ratio still requires a ratio of intractable normalising constants that the authors bypass by Møller et al. (2006) auxiliary variable trick. While the authors mention the alternative pseudo-marginal approach of Andrieu and Roberts (2009), they do not try to implement it, although this would be straightforward here: since the normalising constants are the probabilities of accepting a downward and an upward move, respectively. Those can easily be evaluated at a cost similar to the use of the auxiliary variables. That is,
– generate a few moves from the current value and record the proportion p of accepted downward moves;
– generate a few moves from the final proposed value and record the proportion q of accepted downward moves;
and replace the ratio of intractable normalising constants with p/q. It is not even clear that one needs those extra moves since the algorithm requires an acceptance in the downward and upward moves, hence generate Geometric variates associated with those probabilities p and q, variates that can be used for estimating them. From a theoretical perspective, I also wonder if forcing the downward and upward moves truly leads to an improved convergence speed. Considering the case when the random walk is poorly calibrated for either the downward or upward move, the number of failed attempts before an acceptance may get beyond the reasonable.
If we complement the target π(x) with the conditional P(w|x), the acceptance probability would then involve
{π(x’) P(w’|x’) / π(x) P(w|x)} / {w’ P(w’|x’) / w P(w|x)} = {π(x’) / π(x)} {w/w’}
so it seems the pseudo-marginal (or auxiliary variable) argument also extends to the proposal. Here is a short experiment that shows no discrepancy between target and histogram:
nozero=1e-300 #love-hate move move<-function(x){ bacwa=1;prop1=prop2=rnorm(1,x,2) while (runif(1)>{pi(x)+nozero}/{pi(prop1)+nozero}){ prop1=rnorm(1,x,2);bacwa=bacwa+1} while (runif(1)>{pi(prop2)+nozero}/{pi(prop1)+nozero}) prop2=rnorm(1,prop1,2) y=x if (runif(1)<pi(prop2)*bacwa/pi(x)/fowa){ y=prop2;assign("fowa",bacwa)} return(y)} #arbitrary bimodal target pi<-function(x){.25*dnorm(x)+.75*dnorm(x,mean=5)} #running the chain T=1e5 x=5*rnorm(1);luv8=rep(x,T) fowa=1;prop1=rnorm(1,x,2) #initial estimate while (runif(1)>{pi(x)+nozero}/{pi(prop1)+nozero}){ fowa=fowa+1;prop1=rnorm(1,x,2)} for (t in 2:T) luv8[t]=move(luv8[t-1])
Filed under: Books, pictures, R, Statistics, Travel Tagged: auxiliary variable, doubly intractable problems, Metropolis-Hastings algorithm, Monte Carlo Statistical Methods, multimodality, normalising constant, parallel tempering, pseudo-marginal MCMC, The night of the hunter, unbiased estimation
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.