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)
LDTransformData (AccumulateData Object)
solution 1.808
comb_bound_low 1.808
comb_bound_high 1.809
comb_flags 1
n_total 2^(13)
n 2^(13)
time_integrate 0.017
CubQMCLatticeG (StoppingCriterion Object)
abs_tol 0.001
rel_tol 0
n_init 2^(10)
n_max 2^(35)
CustomFun (Integrand Object)
Gaussian (TrueMeasure Object)
mean 0
covariance 2^(-1)
decomp_type PCA
Lattice (DiscreteDistribution Object)
d 2^(1)
dvec [0 1]
randomize 1
order natural
entropy 273562359450377681412227949180408652150
spawn_key ()
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
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.