rj-MCMC


Summary

This library provides routines for running Reversible Jump Monte-Carlo Markov chains for 1-D and 2-D spatial regression problems (i.e. given noisy x, y data, construct the underlying signal). It also allows generalization to any spatial 1-D and 2-D problem through the users inclusion of a forward model. The routines here are used as the basis of problem specific applications such as Receiver function inversion (code rj-RF) and 2-D travel time tomography (code rj-Tomo). For regression problems the method is also known as Bayesian Partition Modeling (Dension et al. 2002, Gallagher et al. 2011).

The regression problem is treated as one of Bayesian inference and sampled with transdimensional Markov chain. The parametrization of the unknown curve is variable and forms part of the inference process. The x-axis is divided into a set of partitions within each the data is fit with a polynomial with variable order (see Figure). With the package one can re-construct signals from noisy (x,y) data with the number of partitions, the location of the partition boundaries, the order of the polynomial in each partition, and the level of noise on the data (y values) all treated as unknown variables.

Figure of rj-MCMC parametrization for regression problems. Here a zeroth order polynomial is inside 6 partition model (blue curve). The blue nodes indicate that the curve height, width and number of partitions are all variable.

The output of the algorithm is an ensemble of curves with varying numbers of partitions and (optionally) polynomial orders. By taking the mean, median or maximum of the PDF at each x point, an estimated ‘mean’ curve can be constructed from the ensemble of solutions. An example is given in the figure below.


Figure of the mean curve (green) in a 1-D regression problem produced by the rj-MCMC package. The true curve is the piecewise gray line. Red dots are the original data. A histogram of the number of partitions of solutions in the ensemble is shown (right). The true number of partitions is 9 (red line) and algorithm has managed to detect this (From Sambridge et al. 2013)

There are 8 different callable routines in the library, each callable from Fortran90, C or python. Each routine handles a specific problem setup. The Figure below describes the range of generic problems that can be dealt with.


Figure of library routines in rj-MCMC package


Movie of output ensemble for a 2D regression problem


Download

The code can be downloaded here. Enquires should be directed to the current maintainer Julian.Byrne@anu.edu.au. You will need to register with iEarth prior to download.

To unpack the contents of this file, type something like:

    tar xzvf rjmcmc-<version>.tar.gz

The README file provides information on installation, documentation and examples.


References

    Denison, D. G. T., Holmes, C.C., Mallick, B. K., Smith, A. F. M., 2002. Bayesian Methods for Nonlinear Classification and Regression, Wiley, Chichester.

    Gallagher, K., Bodin, T., Sambridge, M., Weiss, D., Kylander, M., and Large, D., 2011. Inference of abrupt changes in noisy geochemical records using Bayesian transdimensional changepoint models, Earth and Planet. Sci. Lett., 311, 182-194, doi:10.1016/j.epsl.2011.09.015

    Sambridge, M., Bodin, T., Gallagher, K., Tkalčić, H., 2013. Transdimensional inference in the geosciences, Phil. Trans. R. Soc. A, 371, 20110547, doi:10.1098/rsta.2011.0547.