Random dive MH

A new Metropolis-Hastings algorithm that I would call “universal” was posted by Somak Dutta yesterday on arXiv. Multiplicative random walk Metropolis-Hastings on the real line contains a different Metropolis-Hastings algorithm called the random dive. The proposed new value x’ given the current value x is defined by

when is a random variable on . Thus, at each iteration, the current value is either shrunk or expanded by a random amount. When I read the paper in the Paris metro, while appreciating that this algorithm could be geometrically ergodic as opposed to the classical random walk algorithm, I was not convinced by the practical impact of the method, as I thought that the special role of zero in the proposal was not always reflected in the target. Especially when considering that the proposal is parameter-free. However, after running the following R program on a target centred away from zero, I found the proposal quite efficient.

f=function(x){.5*dnorm(x,mean=14)+.5*dnorm(x,mean=35)}
Nsim=10^5
x=rep(5,Nsim)
for (t in 2:Nsim){
coef=runif(1,min=-1)^sample(c(-1,1),1)
prop=x[t-1]*coef
prob=abs(coef)*f(prop)/f(x[t-1])
x[t]=x[t-1]
if (runif(1)<prob) x[t]=prop
}
hist(x,pro=T,nclass=113,col=”wheat2″)
curve(f,add=T,n=1001,col=”sienna”,lwd=2)

Obviously, it is difficult to believe that this extension will keep working similarly well when the dimension increases but this is an interesting way of creating a heavy tail proposal.

A repulsive random walk

Matt Asher posted an R experiment on R-bloggers yesterday simulating the random walk

which has the property of avoiding zero by quickly switching to a large value as soon as is small. He was then wondering about the “convergence” of the random walk given that it moves very little once is large enough. The values he found for various horizons t seemed to indicate a stable regime.

I reran the same experiment as Matt in a Monte Carlo perspective, using the R program

resu=matrix(0,ncol=100,nrow=25)
sampl=rnorm(100)
for (i in 1:25){
  for (t in 2^(i-1):2^i) sampl=sampl+rnorm(100)/sampl
     resu[i,]=sampl
     }
boxplot(as.data.frame(t(abs(resu))),name=as.character(1:25),col="wheat3")

The outcome of this R code plotted above shows that the range and the average of the 100 replications is increasing with t. This behaviour indicates a transient behaviour of the Markov chain, which almost surely goes to infinity and never comes back (because at infinity the variance is zero). Another indication for transience is shown by the fact that comes back to the interval (-1,1) with probability , a probability which goes to zero with . As suggested to me by Randal Douc, this transience can be established rigorously by considering

which is thus bounded from below by a null recurrent process, which almost surely goes to infinity. Therefore the above Markov chain cannot have a stationary distribution or even a stationary measure: it almost surely goes to (plus or minus) infinity.

Random mowing

My brother-in-law Christophe has a mowing robot (called Mara!) that is completely autonomous. It starts from its plug and goes along straight lines until it hits an obstacle or the boundaries defined by an underground cable. Then it turns at random and starts again. The concept is quite interesting but I am surprised that it operates on a purely local random walk principle, without any learning about the topology of the terrain. This means that zones with more obstacles are visited less often. (To compensate for that, the mower starts moving in circles each time it hits denser regions of grass.) And that the mower spends a lot of time looking for its plug. An algorithmic determination of the terrain would be possible, obviously, using for instance a Metropolis-Hastings scheme.

t-walk on the banana side

Following my remarks on the t-walk algorithm in the recent A General Purpose Sampling Algorithm for Continuous Distributions, published by Christen and Fox in Bayesian Analysis that acts like a general purpose MCMC algorithm, Darren Wraith tested it on the generic (10 dimension) banana target we used in the cosmology paper. Here is an output from his comparison R program:

The use of the t-walk algorithm (left, with the same number of particles) in this very special case thus produces a wider variability  on the estimated means than the use of adaptive MCMC (center) and our (tuned) PMC algorithm (right).

t-walk on the wild side

When I read in the abstract of the recent A General Purpose Sampling Algorithm for Continuous Distributions, published by Christen and Fox in Bayesian Analysis that

We develop a new general purpose MCMC sampler for arbitrary continuous distributions that requires no tuning.

I am slightly bemused. The proposal of the authors is certainly interesting and widely applicable but to cover arbitrary distributions in arbitrary dimensions with no tuning and great performances sounds too much like marketing on steroids! The 101 Theorem in MCMC methods is that, no matter how good your sampler is, there exists an exotic distribution out there whose only purpose is to make it crash!

The algorithm in A General Purpose Sampling Algorithm for Continuous Distributions is based on two dual and coupled chains which are used towards a double target . Given that only one of the two chains moves at each iteration, according to a random walk, there is a calibration parameter that definitely influences the performances of the method, if not the acceptance probability. This multiple chain approach is reminding me  both of coupled schemes developed by Gareth Roberts in the late 1990′s, along with Laird Breyer, in the wake of the perfect sampling “revolution” and of delayed rejection sampling, as proposed by Antonietta Mira in those years as well. However, there is no particular result in the paper showing an improvement in convergence time over more traditional samplers. (In fact, the random walk nature of the algorithm strongly suggests a lack of uniform ergodicity.) The paper only offers a comparison with an older optimal scaled random walk proposal of Roberts and Rosenthal (Statistical Science, 2001). Rather than with the more recent and effective adaptive Metropolis-Hastings algorithm developed by the same authors.

Since the authors developed a complete set of computer packages, including one in R, I figure people will start to test the method to check for possible improvement over the existing solutions. If the t-walk is indeed superior sui generis, we should hear more about it in the near future…