monte-library is a set of Monte Carlo methods in Python. The package is written to be flexible, clear to understand and encompass variety of Monte Carlo methods.
- Free software: MIT license
- Documentation: https://monte-library.readthedocs.io/en/latest/
Preferred method to install monte-library is through Python's package installer pip. To install monte-library, run this command in your terminal
pip install monte-library
Alternatively, you can install the package directly from GitHub:
git clone -b development https://github.com/draktr/monte-library.git
cd monte-library
python setup.py install
- saving samples and log probability values as
.csvfile - posterior mean, standard deviation and quantiles
- diagnostic tools: effective sample size (ESS), autocorrelation plots, ergodic mean plots, acceptance rate, k-fold histograms, Gelman-Rubin statistic (R-hat), Multivariate Gelman-Rubin Statistic (R-hat), Geweke diagnostic, Heidelberger-Welch diagnostic, Raftery-Lewis diagnostic, Markov chain stationarity tests
- visualizations: histograms, kernel density estimation plots, traceplots
- multidimensional Monte Carlo integration
- multidimensional rejection sampling
- multidimensional importance sampling
- symmetric proposal Metropolis algorithm
- Metropolis-Hastings algorithm
- Gibbs sampler
- vanilla Hamiltonian Monte Carlo
- FLEXIBLE - the package allows users to use various existing Monte Carlo methods for their needs without needing to write the whole algorithm. At the same time,
monte-libraryallows users to specify their own hyperparameters, posterior and proposal distributions as needed. FurthermoreBaseSamplerclass can be used as parent class for any proprietary Monte Carlo algorithm thus utilizing its features for visualizations, posterior analysis and convergence checks. - SIMPLE AND CLEAR CODE BASE - code was intentionally kept simple to be understandable to those with limited exposure to Statistical Computing.
monte-libraryis a great tool to supplement learning as it generally matches mathematical formulations of algorithms and simple syntax helps focus on the algorithm itself. - COMPREHENSIVE - includes Monte Calor methods for various applications. Bayesian modelling methods include both classical methods (Metropolis algorithm) as well as more advanced methods such as Hamiltonian Monte Carlo.
Package contains variety of Monte Carlo methods that can be applied to problems ranging from integration to modelling. Importantly, code is both simple and generalized as to match the respective mathematical formulations of algorithms. As such it can be a great supplement when learning these topics. Finally, the package is flexible and BaseSampler class can be used as a parent class to any user-defined sampler. Furthermore, it is easy to modify existing algorithms with proprietary improvements.
The following example is a simple Monte Carlo implementation to solve the following integral:
from monte import integrator
def integrand(args):
return args[0] ** 2 + args[1] ** 3
result = integrator(integrand, lower_bounds=[-3, -3], upper_bounds=[3, 3], n=10000000)
result
Example 2 is using Metropolis algorithm (with Gaussian proposal) to estimate parameters of a multivariate linear regression.
import numpy as np
from scipy import stats
from monte import GaussianMetropolis
# First, we generate some data
true_theta = np.array([5, 10, 2, 2, 4])
n = 1000
x = np.zeros((n, 4))
x[:, 0] = np.repeat(1, n)
x[:, 1:4] = stats.norm(loc=0, scale=1).rvs(size=(n, 3))
mu = np.matmul(x, true_theta[0:-1])
y = stats.norm(loc=mu, scale=true_theta[-1]).rvs(size=n)
# Define the posterior
def posterior(theta, x, y):
beta_prior = stats.multivariate_normal(
mean=np.repeat(0, len(theta[0:-1])),
cov=np.diag(np.repeat(30, len(theta[0:-1]))),
).logpdf(theta[0:-1])
sd_prior = stats.uniform(loc=0, scale=30).logpdf(theta[-1])
mu = np.matmul(x, theta[0:-1])
likelihood = np.sum(stats.norm(loc=mu, scale=theta[-1]).logpdf(y))
return beta_prior + sd_prior + likelihood
# Lastly, we sample
gaussian_sampler = GaussianMetropolis(posterior)
gaussian_sampler.sample(
iter=10000,
warmup=5000,
theta=np.array([0, 0, 0, 0, 1]),
step_size=1,
lag=1,
x=x,
y=y,
)
Using methods from the BaseSampler class we can perform posterior analytics. These are some of the analytics methods:
# Checking parameter estimates and their credible intervals
gaussian_sampler.mean()
gaussian_sampler.credible_interval()
# Checking Metropolis acceptance rate
gaussian_sampler.acceptance_rate()
# Plotting KDE plot with histogram
gaussian_sampler.parameter_kde()
# Plotting traceplots and ergodic means, and calculating effective sample sizes as convergence diagnostics
gaussian_sampler.traceplots()
gaussian_sampler.plot_ergodic_mean()In the following example we use Hamiltonian Monte Carlo (HMC) algorithm to sample from a distribution. Note that this is a toy example, and HMC is more appropriate to be used for higher-dimensional model parameter estimation. Also note that analytical gradient is not necessary.
import numpy as np
from monte import HamiltonianMC
# Defining the distribution that we are going to sample from...
def posterior(theta):
return -0.5 * np.sum(theta**2)
# ... and its gradient
def posterior_gradient(theta):
return -theta
# Sampling
sampler = HamiltonianMC(posterior, posterior_gradient)
sampler.sample(
iter=10000,
warmup=10,
theta=np.array([8.0, -3.0]),
epsilon=0.01,
l=10,
metric=None,
lag=1,
)There are more sophisticated and computationally efficient implementations of Monte Carlo methods for off-the-shelf solutions
- ArviZ - independent library for exploratory analysis of Bayesian models
- vegas - uses improved version of the adaptive Monte Carlo vegas algorithm
- OpenBUGS - open source implementation of BUGS language utilizing Gibbs sampler
- JAGS - cross-platform and more extensible implementation of BUGS language
- WinBUGS - software for Bayesian analysis utilizing Gibbs sampler (available, but discontinued in favour of OpenBUGS)
- Stan - state-of-the-art probabilistic language implementing advanced version of No-U-Turn Sampler
- PyMC - supports HMC and Metropolis-Hastings algorithms, as well as Sequential Monte Carlo methods
- Easy to be understood and used by non-mathematicians
- Potential to be used as pedagogical tool
- Easy to modify algorithms with proprietary improvements
- Flexibility and simplicity over computational efficiency
- Tested
- Dedicated documentation
- Formatting deferred to Black
Feel free to reach out through Issues forum if you wish to add features or help in any other way. If there are any issues, bugs or improvement recommendations, please raise them in the forum. Especially reach out if you want to contribute with any of the possible features listed below.
- Sampling trace animation
- ECDF plot
- Forrest plot of parameter estimates with credible intervals
- Monte Carlo Error
- Support for multiple chains
- Slice sampling
- Annealed importance sampling
- Component-wise Metropolis
- Independent Metropolis
- Wang-Landau algorithm
- Monte Carlo tree search
- Direct Sampling Monte Carlo
- Monte Carlo statistical distribution test
The following is the non-exhaustive list of useful sources for learning more about Monte Carlo methods. Some of the code in monte has been written based on mathematical formulae from some of these sources.
[1] Ntzoufras, I. (2009). Bayesian Modelling Using WinBUGS. Wiley.
[2] Metropolis, N., & Ulam, S. (1949). The Monte Carlo Method. Journal of the American Statistical Association, 44(247), 335–341. https://doi.org/10.1080/01621459.1949.10483310
[3] Metropolis, N., Rosenbluth, A. W., Rosenbluth, M. N., Teller, A. H., & Teller, E. (2004). Equation of State Calculations by Fast Computing Machines. The Journal of Chemical Physics, 21(6), 1087. https://doi.org/10.1063/1.1699114
[4] Hastings, W. K. (1970). Monte Carlo sampling methods using Markov chains and their applications. Biometrika, 57(1), 97–109. https://doi.org/10.1093/BIOMET/57.1.97
[5] Hartig, F. (n.d.). A simple Metropolis-Hastings MCMC in R | theoretical ecology. Retrieved February 15, 2023, from https://theoreticalecology.wordpress.com/2010/09/17/metropolis-hastings-mcmc-in-r/
[6] Dirty Quant @YouTube. (n.d.). The Metropolis-Hastings Algorithm (MCMC in Python) - YouTube. Retrieved February 15, 2023, from https://www.youtube.com/watch?v=MxI78mpq_44
[7] TWEAG Software Innovation Lab. (n.d.). Markov chain Monte Carlo (MCMC) Sampling, Part 1: The Basics - Tweag. Retrieved February 15, 2023, from https://www.tweag.io/blog/2019-10-25-mcmc-intro1/
[8] Urbanevych, V. (n.d.). VU | Bayesian linear regression and Metropolis-Hastings sampler. Retrieved February 15, 2023, from https://vitaliiur.github.io/blog/2021/linreg/
[9] Geman, S., & Geman, D. (1984). Stochastic Relaxation, Gibbs Distributions, and the Bayesian Restoration of Images. IEEE Transactions on Pattern Analysis and Machine Intelligence, PAMI-6(6), 721–741. https://doi.org/10.1109/TPAMI.1984.4767596
[10] Campbell, K. R. (n.d.). Gibbs sampling for Bayesian linear regression in Python | Kieran R Campbell - blog. Retrieved February 15, 2023, from https://kieranrcampbell.github.io/blog/2016/05/15/gibbs-sampling-bayesian-linear-regression.html
[11] Betancourt, M. (2017). A Conceptual Introduction to Hamiltonian Monte Carlo. https://doi.org/10.48550/arxiv.1701.02434
[12] Neal, R. M. (2012). MCMC using Hamiltonian dynamics. Handbook of Markov Chain Monte Carlo, 1–592. https://doi.org/10.1201/b10905
[13] Stan. (n.d.). 15.1 Hamiltonian Monte Carlo | Stan Reference Manual. Retrieved February 15, 2023, from https://mc-stan.org/docs/reference-manual/hamiltonian-monte-carlo.html
[14] Clark, M. (n.d.). Hamiltonian Monte Carlo | Model Estimation by Example. Retrieved February 15, 2023, from https://m-clark.github.io/models-by-example/hamiltonian-monte-carlo.html
[15] Richard. (n.d.). Markov Chains: Why Walk When You Can Flow? | Elements of Evolutionary Anthropology. Retrieved February 15, 2023, from http://elevanth.org/blog/2017/11/28/build-a-better-markov-chain/