*This work was originally published here, and has been reproduced with the author’s permission.*

We have created QMCPy, an open-source, object-oriented quasi- Monte Carlo (QMC) software framework in Python 3, which contains standardized parent classes for modeling integrands, measure, discrete distribution, and stopping criteria. We hope QMCPy could enable researchers and users to quickly extend and experiment with novel algorithmic components that validate their theories, or to simply leverage or compare proven ones as they develop their applications.

In this tutorial, we introduce QMCPy [1] by an example. QMCPy can be installed with `pip install qmcpy`

or cloned from the QMCSoftware GitHub repository.

Consider the problem of integrating the Keister function [2] with respect to a \(d\)-dimensional Gaussian measure:

\[f(\boldsymbol{x}) = \pi^{d/2} \cos(||\boldsymbol{x}||), \qquad \boldsymbol{x} \in \mathbb{R}^d, \qquad \boldsymbol{X} \sim \mathcal{N}(\boldsymbol{0}_d,\boldsymbol{I}_d/2),\] \[\begin{aligned} \mu &= \mathbb{E}[f(\boldsymbol{X})] := \int_{\mathbb{R}^d} f(\boldsymbol{x}) \, \pi^{-d/2} \exp( – ||\boldsymbol{x}||^2) \, \rm d \boldsymbol{x} \\ &= \int_{[0,1]^d} \pi^{d/2} \cos\left(\sqrt{ \frac{1}{2} \sum_{j=1}^d\Phi^{-1}(x_j)}\right) \, \rm d \boldsymbol{x} \end{aligned}\]

where \(\|\boldsymbol{x}\|\) is the Euclidean norm, \(\boldsymbol{I}_d\) is the \(d\)-dimensional identity matrix, and \(\Phi\) denotes the standard normal cumulative distribution function. When \(d=2\), \(\mu \approx 1.80819\) and we can visualize the Keister function and realizations of the sampling points depending on the tolerance values, \(\varepsilon\), in the following figure.

The Keister function is implemented below with help from NumPy [3] in the following code snippet:

import numpy as np def keister(x): """ x: nxd numpy ndarray with n samples d dimensions returns n-vector of the Kesiter function evaluations """ d = x.shape[1] norm_x = np.sqrt((x**2).sum(1)) k = np.pi**(d/2) * np.cos(norm_x) return k # size n vector

In addition to our Keister integrand and Gaussian true measure, we must select a discrete distribution and a stopping criterion [4]. The stopping criterion determines the number of points at which to evaluate the integrand in order for the mean approximation to be accurate within a user-specified error tolerance, $\varepsilon$. The discrete distribution determines the sites at which the integrand is evaluated.

For this Keister example, we select the lattice sequence as the discrete distribution and corresponding cubature-based stopping criterion [5]. The discrete distribution, true measure, integrand, and stopping criterion are then constructed within the QMCPy framework below.

import qmcpy lattice = qmcpy.Lattice(dimension = 2) gaussian = qmcpy.Gaussian(lattice, mean = 0, covariance = 1/2) cf_keister = qmcpy.CustomFun(gaussian, custom_fun = keister) stopping_criterion = qmcpy.CubQMCLatticeG(cf_keister, abs_tol = 1e-4)

Calling `integrate`

on the `stopping_criterion`

instance returns the numerical solution and a data object. Printing the data object will provide a neat summary of the integration problem. For details of the output fields, refer to the online, searchable QMCPy Documentation at https://qmcpy.readthedocs.io.

solution, data = stopping_criterion.integrate() print(data)

```
Solution: 1.8082
CustomFun (Integrand Object)
Lattice (DiscreteDistribution Object)
dimension 2^(1)
randomize 1
seed None
backend gail
mimics StdUniform
Gaussian (TrueMeasure Object)
mean 0
covariance 2^(-1)
CubQMCLatticeG (StoppingCriterion Object)
abs_tol 1.00e-04
rel_tol 0
n_init 2^(10)
n_max 2^(35)
LDTransformData (AccumulateData Object)
n_total 2^(16)
solution 1.808
r_lag 2^(2)
time_integrate 0.070
```

This guide is not meant to be exhaustive but rather a quick introduction to the QMCPy framework and syntax. In an upcoming blog, we will take a closer look at low-discrepancy sequences such as the lattice sequence from the above example.

#### References

1. Choi, S.-C. T., Hickernell, F. J., McCourt, M., Rathinavel, J., & Sorokin, A. QMCPy: A quasi-Monte Carlo Python Library, https://qmcsoftware.github.io/QMCSoftware/. 2020.

2. Keister, B. D., Multidimensional Quadrature Algorithms. Computers in Physics 10, 119–122 (1996).

3. Oliphant, T., Guide to NumPy, https://ecs.wgtn.ac.nz/foswiki/pub/Support/ManualPagesAndDocumentation/numpybook.pdf (Trelgol Publishing, USA, 2006).

4. Hickernell, F. J., Choi, S.-C. T., Jiang, L. & Jimenez Rugama, L. A. in WileyStatsRef-Statistics Reference Online (eds Davidian, M. et al.) (John Wiley & Sons Ltd., 2018).

5. Jimenez Rugama, L. A. & Hickernell, F. J., Adaptive Multidimensional Integration Based on Rank-1 Lattices in Monte Carlo and Quasi-Monte Carlo Methods: MCQMC, Leuven, Belgium, April 2014 (eds Cools, R. & Nuyens, D.) 163. arXiv:1411.1966 (Springer-Verlag, Berlin, 2016), 407–422.

Use SigOpt free. Sign up today.