A QMCPy Quick Start¶
In this tutorial, we introduce QMCPy  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  with respect to a -dimensional Gaussian measure:
where is the Euclidean norm, is the -dimensional identity matrix, and denotes the standard normal cumulative distribution function. When , and we can visualize the Keister function and realizations of the sampling points depending on the tolerance values, , in the following figure:
The Keister function is implemented below with help from NumPy  in the following code snippet:
import numpy as np def keister(x): """ x: nxd numpy ndarray n samples d dimensions returns n-vector of the Kesiter function evaluated at the n input samples """ d = x.shape 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 . 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, . 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 . The discrete distribution, true measure, integrand, and stopping criterion are then constructed within the QMCPy framework below.
import qmcpy d = 2 discrete_distrib = qmcpy.Lattice(dimension = d) true_measure = qmcpy.Gaussian(discrete_distrib, mean = 0, covariance = 1/2) integrand = qmcpy.CustomFun(true_measure,keister) stopping_criterion = qmcpy.CubQMCLatticeG(integrand = integrand, abs_tol = 1e-3)
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) d 2^(1) randomize 1 order natural seed 476604 mimics StdUniform Gaussian (TrueMeasure Object) mean 0 covariance 2^(-1) decomp_type pca CubQMCLatticeG (StoppingCriterion Object) abs_tol 0.001 rel_tol 0 n_init 2^(10) n_max 2^(35) LDTransformData (AccumulateData Object) n_total 2^(13) solution 1.808 error_bound 5.93e-04 time_integrate 0.016
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.
Choi, S.-C. T., Hickernell, F., McCourt, M., Rathinavel J., & Sorokin, A. QMCPy: A quasi-Monte Carlo Python Library. https://qmcsoftware.github.io/QMCSoftware/. 2020.
Keister, B. D. Multidimensional Quadrature Algorithms. Computers in Physics 10, 119–122 (1996).
Oliphant, T., Guide to NumPy https://ecs.wgtn.ac.nz/foswiki/pub/Support/ManualPagesAndDocumentation/numpybook.pdf (Trelgol Publishing USA, 2006).
Hickernell, F., 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).
Jimenez Rugama, L. A. & Hickernell, F. Adaptive Multidimensional Inte-gration 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.