Monte Library

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/

Installation

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

Features

Base module

saving samples and log probability values as .csv file
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

General Monte Carlo Methods

multidimensional Monte Carlo integration
multidimensional rejection sampling
multidimensional importance sampling

Markov Chain Monte Carlo Modelling Methods

symmetric proposal Metropolis algorithm
Metropolis-Hastings algorithm
Gibbs sampler
vanilla Hamiltonian Monte Carlo

Advantages

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-library allows users to specify their own hyperparameters, posterior and proposal distributions as needed. Furthermore BaseSampler class 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-library is 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.

Usage

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.

Example 1: Monte Carlo Integration

The following example is a simple Monte Carlo implementation to solve the following integral: $$\int_{-3}^{3} \int_{-3}^{3} x^2 + y^3 dxdy$$

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: Bayesian Linear Regression with Metropolis Algorithm

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()

Example 3: Sampling from a Multivariate Distribution using Hamiltonian Monte Carlo

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,
    )

Alternatives and Complements

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

Project Principles

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

Future Development

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.

Possible Future Features

In `BaseSampler` Class

Sampling trace animation
ECDF plot
Forrest plot of parameter estimates with credible intervals
Monte Carlo Error
Support for multiple chains

Monte Carlo Methods

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

[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/

Gibbs Sampler

[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

Hamiltonian Monte Carlo

[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/

draktr/monte-library