Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.
target=function(x,a=1,b=1){dbeta(x,a,b)} T=1e4 #MCMC runs bsch=50 #batsch size mkv=rep(runif(1),T) #MCMC chain eps=rep(1/runif(1,10,100),T) #adapted proposal scale ace=rep(0,T) #acceptance rate for (t in 2:T){ batsch=mkv[t-1] for (d in 1:bsch){ # Beta proposal prop=rbeta(1,eps[t-1]*batsch,eps[t-1]*(1-batsch)) mh=target(prop)*dbeta(batsch,eps[t-1]* prop,eps[t-1]*(1-prop))/ dbeta(prop,eps[t-1]*batsch,eps[t-1]*(1-batsch))/ target(batsch) #Metropolis-Hastings ratio if (runif(1)<mh){ace[t-1]=ace[t-1]+1;batsch=prop}} ace[t-1]=ace[t-1]/bsch mkv[t]=batsch #subsampled Markov chain # update of epsilon if (ace[t-1]<.456){eps[t]=eps[t-1]*1.01 }else{eps[t]=eps[t-1]/1.01} } #end of the MCMC loop
Indeed, when using a Be(a,b) target with either a<1 or b<1, the above MCMC algorithm will eventually but inevitably run into difficulties, because the Beta generator ends up producing either 0 or 1, thanks to rounding effects, hence an NA value for the Metropolis-Hastings ratio. For a Be(a,b) target with both coefficients above one, there was no convergence issue and no difficulty with the behaviour of ε. As shown by the picture above when the target is uniform. Or the one below when it is a Be(4,8).
Filed under: Books, Kids, R, Statistics, University life Tagged: adaptive MCMC methods, beta distribution, Dirichlet prior, Gaussian random walk, Markov chains
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.