austerity in MCMC land (#2)

After 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]}
 }

interesting puzzle

In addition to its weekly mathematics puzzles, Le Monde is now publishing a series of vulgarisation books on mathematics, under the patronage of Cédric Villani. Jean-Michel Marin brought me two from the series, one on the golden number and one on Pythagoras’ theorem. (This is actually a translation of a series published by El Pais last year.) Those books are a bit stretched given the topic, even though I enjoyed the golden number (the second one having a lot of redundancy with the first one.) However, I came upon an interesting question, namely about the maximum size of a cube that could fit through a tunnel drilled through the unit cube. Sadly, I could not find an answer to this problem on the web, even though the book mentions a solution with a side larger than one…

Le Monde puzzle [#817]

The weekly Le Monde puzzle is (again) a permutation problem that can be rephrased as follows:

Find

where denotes the set of permutations on {0,…,10} and is defined modulo 11 [to turn {0,...,10} into a torus]. Same question for

and for

This is rather straightforward to code if one adopts a brute-force approach::

perminmax=function(T=10^3){
  mins=sums=rep(500,3)
  permin=matrix(0:10,ncol=11,nrow=3,byrow=TRUE)

  for (t in 1:T){
    perms=sample(0:10)
    adde=perms+perms[c(11,1:10)]
    sums[1]=max(adde)
    adde=adde+perms[c(10,11,1:9)]
    sums[2]=max(adde)
    adde=adde+perms[c(9:11,1:8)]+perms[c(8:11,1:7)]
    sums[3]=max(adde)
    for (j in 1:3)
     if (sums[j]<mins[j]){
       mins[j]=sums[j];permin[j,]=perms}
    }

  print(mins)
  print(permin)
  }

(where I imposed the first term to be zero because of the invariance by permutation), getting the solutions

> perminmax(10^5)
[1] 11 17 28
     [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11]
[1,]    0   10    1    6    5    4    7    3    8     2     9
[2,]    0   10    4    3    5    1    9    6    2     8     7
[3,]    0    2    9    6    7    3    1    4   10     8     5

for 2, 3, and 5 terms.  (Since 10! is a mere 3.6 million, I checked with an exhaustive search, using the function permutation from the gtools package.)

MCMSki IV, Jan. 6-8, 2014, Chamonix (news #5)

More exciting news about MCMSki IV!

First thing first, the 16 contributed sessions are now all-set, having gotten the stamp of approval from the scientific committee! Thanks to everyone who submitted a session proposal. (There were so many proposals that we alas had to reject some, as well as every single talk proposal… Sorry people: we hope to hear about your research advances via your posters!) See the MCMSki IV website for the whole list. Apart from the plenary lectures, and the round table on software held on the second evening, there will be three parallel sessions on the remaining three slots for each day of the conference, which means 25 sessions total!

Second, the “call for posters” is open, simply meaning that anyone wishing to present a poster at MCMSki IV on Monday evening (or Tuesday night if we cannot accommodate all posters within a single evening!) is welcome to do so! This will take place in the conference centre as well (with an open bar to keep up with traditions) To this effect, if you intend to present a poster, (a) tick the box in the registration form and (b) …wait for further instructions on the MCMSki IV website about sending your abstract as we are trying to find an easy way to store and publish posters there. Simple as AB(C)!

Last, the registration page is now open! So fell free to register at your earliest convenience. The deadline for early bird registration is October 15, 2013 however hotel rooms are likely to vanish much earlier than that, leaving you on your own to find accommodation in Chamonix (not such a terrible task, actually!)

BayesComp homepage

Today, the BayesComp section of ISBA launched its website. It is organised as a wiki and members of the section are strongly incited to take part into the construction of the website. To quote from Peter Green’s introduction:

This new Wikidot site aims to be a community-edited resource on all aspects of Bayesian computation, available for all to read; here ‘community’ means members of the section – we hope that interested members will help us create pages of information and advice to help disseminate new research ideas in an accessible way, and promote good practice. Members can also submit links to a directory of papers, slides, videos and software, and to a diary of upcoming events such as conferences and workshops, all through easy-to-use forms.

Please visit BayesComp, read the Quick Start Guide there, actively contribute to the site’s content, and let us have feedback using the blog facility provided. Once you have editing credentials on the site, you do not need permission to make edits, just do it! Section officers may tidy things up later, if no one else does, but we won’t delete anything unless it is offensive or plainly wrong.

Thanks to the inputs of Peter Green and of Nicolas Chopin, this could be a wonderful exchange tool for the community, but only if this community strives to keep it alive!