trying to speed up Metropolis… and failing!

June 12, 2014
By

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

A while ago (but still after Iceland since I used the thorn rune as a math symbol!), I wrote the following post draft as a memo. Now that Marco Banterle, Clara Grazian and myself have completed our delayed acceptance paper, it may be of interest to some readers to see how a first attempt proved fruitless.

In the past days, I tried to speed up my student Clara’s code, reducing the number of prior evaluations in a Metropolis-Hastings algorithm by trying to reject proposals with low likelihoods before computing the corresponding prior. (This is one of those rare problems where the prior is the culprit.) My first idea was to start from the case when the new likelihood þ(θ’) was lower than the previous likelihood þ(θ)

þ(θ’) < þ(θ)

as the indicator

ℑ{u≤þ(θ’)/þ(θ)}

is an unbiased estimator of the ratio þ(θ’)/þ(θ) . When u is generated from a uniform U(0,1). (All u’s in this page will be uniform U(0,1), unless otherwise indicated.) Does this mean low values of þ(θ’) can be rejected prior to computing π(θ’)? Unfortunately no, since

ℑ{u≤þ(θ’)/þ(θ)} π(θ’)/ π(θ) ∧ 1

is not an unbiased estimator of

{þ(θ’)/þ(θ)} {π(θ’)/ π(θ)} ∧ 1

Then I considered

{þ(θ’)/þ(θ)} {π(θ’)/ π(θ)} ∧ þ(θ)/þ(θ’)}

which can be associated with the unbiased estimator

ℑ{u≤þ(θ’)/þ(θ)} {π(θ’)/ π(θ) ∧ þ(θ)/þ(θ’)}

Once more unfortunately, since

{π(θ’)/ π(θ) ∧ þ(θ)/þ(θ’)}

is not easily bounded, unless π(θ’)/π(θ) itself is bounded (a self-defeating condition!), generating an unbiased estimator of this second term is not obvious.

delayedAnother attempt was to consider an exact approximation à la Andrieu and Roberts (2009, Annals of Stat.), using an unbiased estimator of þ(θ’), namely

ℑ{þ(θ)u≤þ(θ’)} þ(θ) ℑ{þ(θ)>þ(θ’)} + þ(θ’) ℑ{þ(θ)<þ(θ’)}

and to plug this unbiased estimator in the Metropolis-Hastings acceptance ratio,

π(θ’) [ℑ{þ(θ)u’≤þ(θ’)}þ(θ)ℑ{þ(θ)>þ(θ’)} + þ(θ’)ℑ{þ(θ)<þ(θ’)}] / π(θ) [ℑ{þ(θ’)u≤þ(θ)}þ(θ’)ℑ{þ(θ’)>þ(θ)} + þ(θ)ℑ{þ(θ)>þ(θ’)}] ∧ 1

which leads to

π(θ’) [ℑ{þ(θ)u’≤þ(θ’)}þ(θ)] / π(θ) þ(θ) ∧ 1 = π(θ’) ℑ{þ(θ)u’≤þ(θ’)} / π(θ) ∧ 1

when þ(θ)>þ(θ’) and to

π(θ’) þ(θ’) / π(θ) ℑ{þ(θ’)u≤þ(θ)}þ(θ’) ∧ 1 = π(θ’) / π(θ) ℑ{þ(θ’)u≤þ(θ)} ∧ 1

otherwise. But this is not correct either. As shown by the above histogram versus target.

#data  x~B(N,p)
#posterior p|x~Be(x+7.5,N+.5-x)

T=10^4
N=100;y=32
a=7.5;b=.5

chk=NULL
p=rep(.2,T)
for (t in 2:T){

  pp=p[t-1]+runif(1,-.1,.1) #random walk
  p[t]=p[t-1] #by default rejection

  if ((pp>0)&(pp<1)){

   ratio=dbinom(y,N,pp)/dbinom(y,N,p[t-1])
   if (ratio<1){

     if ((runif(1)1)||(runif(1)<
         dbeta(pp,a,b)/dbeta(p[t-1],a,b)))
                p[t]=pp}
    }}

Filed under: R, Statistics, University life Tagged: big data, delayed acceptance, pseudo-marginal MCMC, unbiased estimation

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

R-bloggers.com offers daily e-mail updates about R news and tutorials on topics such as: Data science, Big Data, R jobs, 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.

Sponsors

Never miss an update!
Subscribe to R-bloggers to receive
e-mails with the latest R posts.
(You will not see this message again.)

Click here to close (This popup will not appear again)