# A QMCPy Quick Start¶

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 -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:

Keister Function

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
n samples
d dimensions

returns n-vector of the Kesiter function
evaluated at the n input samples
"""
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, . 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
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.

## References¶

1. 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.

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