austerity in MCMC land (#2)

April 29, 2013
By

(This article was first published on Xi'an's Og » R, and kindly contributed to R-bloggers)

mcmc run with median instead of meanAfter reading the arXiv paper by Korattikara, Chen and Welling, I wondered about the expression of the acceptance step of the Metropolis-Hastings algorithm as a mean of log-likelihoods over the sample. More specifically the long sleepless nights at the hospital led me to ponder the rather silly question of the impact of replacing mean by median. I thus tried running a Metropolis-Hastings algorithm with the substitute and it (of course!) let to a nonsensical answer, as shown by the above graph. The true posterior is the one for a normal model and the histogram indicates a lack of convergence of the Markov chain to this posterior even though it does converge to some posterior. Here is the R code for this tiny experiment:

#data generation
N=100
x=rnorm(N)

#HM steps
T=10^5
theta=rep(0,T)
curlike=dnorm(x,log=TRUE)
for (t in 2:T){

  prop=theta[t-1]+.1*rnorm(1)
  proplike=dnorm(x,mean=prop,log=TRUE)
  u=runif(1)
  bound=log(u)-dnorm(prop,sd=10,log=TRUE)+
         dnorm(theta[t-1],sd=10,log=TRUE)
  if (median(proplike)-median(curlike)>bound/N){
   theta[t]=prop;curlike=proplike
   } else { theta[t]=theta[t-1]}
 }

Filed under: R, Statistics Tagged: arXiv, mean vs. median, Metropolis-Hastings, R

To leave a comment for the author, please follow the link and comment on his blog: Xi'an's Og » R.

R-bloggers.com offers daily e-mail updates about R news and tutorials on topics such as: visualization (ggplot2, Boxplots, maps, animation), programming (RStudio, Sweave, LaTeX, SQL, Eclipse, git, hadoop, Web Scraping) statistics (regression, PCA, time series, trading) and more...



If you got this far, why not subscribe for updates from the site? Choose your flavor: e-mail, twitter, RSS, or facebook...

Comments are closed.