Python bindings for Monte Carlo Markov chain (MCMC) diagnostics using the coda package in R
- Install R (https://www.r-project.org/) and rpy2
On Mac with Homebrew
brew tap homebrew/science
brew install r
pip install rpy2
or Ubuntu
sudo apt-get install r-base r-base-dev python-rpy2
- Install the coda package
mkdir $HOME/Rpackages
echo R_LIBS_USER="$HOME/Rpackages" > ~/.Renviron
Rscript -e 'install.packages("coda", repos="http://cran.us.r-project.org")'
- Install py-coda
git clone https://github.com/surhudm/py-coda.git
cd py-coda
python setup.py install
Prepare MCMC file say chainfile.out where every line is the parameter set at
every step of the chain. Prepare another file where each line lists the
parameter name, write this parameter name in latex form, as it will get set as a
label on the plots. For examples, see the files chainfile.ind
and
chainfile.out
(the latter file only contains 100 lines to prevent the
repository from being bulky). Then you could run the diagnostics as follows:
from py_coda import read_pycoda
mcmcobj = read_pycoda("chainfile.out", "chainnew.ind", thin=1)
mcmcobj.geweke()
mcmcobj.get_stats()
mcmcobj.plot_traces()
mcmcobj.plot_autocorr()
mcmcobj.heidelberger_welch()
In your python code you could even use the python bindings for coda directly. For example, if you have a numpy array with shape (Nchain, Nparam) which stores the MCMC where Nchain is the number of iterations of the MCMC and Nparam is the number of parameters, and an array of labels for your parameters, you could simply do
from py_coda import mcmc
mcmcobj = mcmc(labels, data, thin=1)
mcmcobj.geweke()
mcmcobj.get_stats()
mcmcobj.plot_traces()
mcmcobj.plot_autocorr()
mcmcobj.heidelberger_welch()
Please file bug reports or issues if you find any.