QMCPy Documentation

_images/discrete_distribution_overview.png _images/true_measure_overview.png _images/integrand_overview.png _images/stopping_criterion_overview.png

Discrete Distribution Class

_images/discrete_distribution_specific.png

Abstract Discrete Distribution Class

class qmcpy.discrete_distribution._discrete_distribution.DiscreteDistribution(dimension, seed)

Discrete Distribution abstract class. DO NOT INSTANTIATE.

__init__(dimension, seed)
Parameters
  • dimension (int or ndarray) –

    dimension of the generator. If an int is passed in, use sequence dimensions [0,…,dimensions-1]. If a ndarray is passed in, use these dimension indices in the sequence.

    Note that this is not relevent for IID generators.

  • seed (int or numpy.random.SeedSequence) – seed to create random number generator

gen_samples(*args)

ABSTRACT METHOD to generate samples from this discrete distribution.

Parameters

args (tuple) – tuple of positional argument. See implementations for details

Returns

n x d array of samples

Return type

ndarray

pdf(x)

ABSTRACT METHOD to evaluate pdf of distribution the samples mimic at locations of x.

spawn(s=1, dimensions=None)

Spawn new instances of the current discrete distribution but with new seeds and dimensions. Developed for multi-level and multi-replication (Q)MC algorithms.

Parameters
  • s (int) – number of spawn

  • dimensions (ndarray) – length s array of dimension for each spawn. Defaults to current dimension

Returns

list of DiscreteDistribution instances with new seeds and dimensions

Return type

list

Digital Net Base 2

class qmcpy.discrete_distribution.digital_net_b2.digital_net_b2.DigitalNetB2(dimension=1, randomize='LMS_DS', graycode=False, seed=None, generating_matrices='sobol_mat.21201.32.32.msb.npy', d_max=None, t_max=None, m_max=None, msb=None, t_lms=None, _verbose=False)

Quasi-Random digital nets in base 2.

>>> dnb2 = DigitalNetB2(2,seed=7)
>>> dnb2.gen_samples(4)
array([[0.56269008, 0.17377997],
       [0.346653  , 0.65070632],
       [0.82074548, 0.95490574],
       [0.10422261, 0.49458097]])
>>> dnb2.gen_samples(1)
array([[0.56269008, 0.17377997]])
>>> dnb2
DigitalNetB2 (DiscreteDistribution Object)
    d               2^(1)
    dvec            [0 1]
    randomize       LMS_DS
    graycode        0
    entropy         7
    spawn_key       ()
>>> DigitalNetB2(dimension=2,randomize=False,graycode=True).gen_samples(n_min=2,n_max=4)
array([[0.75, 0.25],
       [0.25, 0.75]])
>>> DigitalNetB2(dimension=2,randomize=False,graycode=False).gen_samples(n_min=2,n_max=4)
array([[0.25, 0.75],
       [0.75, 0.25]])

References

[1] Marius Hofert and Christiane Lemieux (2019). qrng: (Randomized) Quasi-Random Number Generators. R package version 0.0-7. https://CRAN.R-project.org/package=qrng.

[2] Faure, Henri, and Christiane Lemieux. “Implementation of Irreducible Sobol’ Sequences in Prime Power Bases.” Mathematics and Computers in Simulation 161 (2019): 13–22. Crossref. Web.

[3] F.Y. Kuo & D. Nuyens. Application of quasi-Monte Carlo methods to elliptic PDEs with random diffusion coefficients - a survey of analysis and implementation, Foundations of Computational Mathematics, 16(6):1631-1696, 2016. springer link: https://link.springer.com/article/10.1007/s10208-016-9329-5 arxiv link: https://arxiv.org/abs/1606.06613

[4] D. Nuyens, The Magic Point Shop of QMC point generators and generating vectors. MATLAB and Python software, 2018. Available from https://people.cs.kuleuven.be/~dirk.nuyens/ https://bitbucket.org/dnuyens/qmc-generators/src/cb0f2fb10fa9c9f2665e41419097781b611daa1e/cpp/digitalseq_b2g.hpp

[5] Paszke, A., Gross, S., Massa, F., Lerer, A., Bradbury, J., Chanan, G., … Chintala, S. (2019). PyTorch: An Imperative Style, High-Performance Deep Learning Library. In H. Wallach, H. Larochelle, A. Beygelzimer, F. d extquotesingle Alch'e-Buc, E. Fox, & R. Garnett (Eds.), Advances in Neural Information Processing Systems 32 (pp. 8024–8035). Curran Associates, Inc. Retrieved from http://papers.neurips.cc/paper/9015-pytorch-an-imperative-style-high-performance-deep-learning-library.pdf

[6] I.M. Sobol’, V.I. Turchaninov, Yu.L. Levitan, B.V. Shukhman: “Quasi-Random Sequence Generators” Keldysh Institute of Applied Mathematics, Russian Acamdey of Sciences, Moscow (1992).

[7] Sobol, Ilya & Asotsky, Danil & Kreinin, Alexander & Kucherenko, Sergei. (2011). Construction and Comparison of High-Dimensional Sobol’ Generators. Wilmott. 2011. 10.1002/wilm.10056.

[8] Paul Bratley and Bennett L. Fox. 1988. Algorithm 659: Implementing Sobol’s quasirandom sequence generator. ACM Trans. Math. Softw. 14, 1 (March 1988), 88–100. DOI:https://doi.org/10.1145/42288.214372

__init__(dimension=1, randomize='LMS_DS', graycode=False, seed=None, generating_matrices='sobol_mat.21201.32.32.msb.npy', d_max=None, t_max=None, m_max=None, msb=None, t_lms=None, _verbose=False)
Parameters
  • dimension (int or ndarray) – dimension of the generator. If an int is passed in, use sequence dimensions [0,…,dimensions-1]. If a ndarray is passed in, use these dimension indices in the sequence.

  • randomize (bool) – apply randomization? True defaults to LMS_DS. Can also explicitly pass in ‘LMS_DS’: linear matrix scramble with digital shift ‘LMS’: linear matrix scramble only ‘DS’: digital shift only

  • graycode (bool) – indicator to use graycode ordering (True) or natural ordering (False)

  • seed (list) – int seed of list of seeds, one for each dimension.

  • generating_matrices (ndarray or str) – generating matricies or path to generating matrices. ndarray should have shape (d_max, m_max) where each int has t_max bits generating_matrices sould be formatted like gen_mat.21201.32.32.msb.npy with name.d_max.t_max.m_max.{msb,lsb}.npy

  • d_max (int) – max dimension

  • t_max (int) – number of bits in each int of each generating matrix. aka: number of rows in a generating matrix with ints expanded into columns

  • m_max (int) – 2^m_max is the number of samples supported. aka: number of columns in a generating matrix with ints expanded into columns

  • msb (bool) – bit storage as ints. e.g. if t_max=3, then 6 is [1 1 0] in MSB (True) and [0 1 1] in LSB (False)

  • t_lms (int) – LMS scrambling matrix will be t_lms x t_max for generating matrix of shape t_max x m_max

  • _verbose (bool) – print randomization details

gen_samples(n=None, n_min=0, n_max=8, warn=True, return_unrandomized=False)

Generate samples

Parameters
  • n (int) – if n is supplied, generate from n_min=0 to n_max=n samples. Otherwise use the n_min and n_max explicitly supplied as the following 2 arguments

  • n_min (int) – Starting index of sequence.

  • n_max (int) – Final index of sequence.

  • return_unrandomized (bool) – return unrandomized samples as well? If True, return randomized_samples,unrandomized_samples. Note that this only applies when randomize includes Digital Shift. Also note that unrandomized samples included linear matrix scrambling if applicable.

Returns

(n_max-n_min) x d (dimension) array of samples

Return type

ndarray

pdf(x)

pdf of a standard uniform

class qmcpy.discrete_distribution.digital_net_b2.digital_net_b2.Sobol(dimension=1, randomize='LMS_DS', graycode=False, seed=None, generating_matrices='sobol_mat.21201.32.32.msb.npy', d_max=None, t_max=None, m_max=None, msb=None, t_lms=None, _verbose=False)

Lattice

class qmcpy.discrete_distribution.lattice.lattice.Lattice(dimension=1, randomize=True, order='natural', seed=None, generating_vector='lattice_vec.3600.20.npy', d_max=None, m_max=None)

Quasi-Random Lattice nets in base 2.

>>> l = Lattice(2,seed=7)
>>> l.gen_samples(4)
array([[0.04386058, 0.58727432],
       [0.54386058, 0.08727432],
       [0.29386058, 0.33727432],
       [0.79386058, 0.83727432]])
>>> l.gen_samples(1)
array([[0.04386058, 0.58727432]])
>>> l
Lattice (DiscreteDistribution Object)
    d               2^(1)
    dvec            [0 1]
    randomize       1
    order           natural
    entropy         7
    spawn_key       ()
>>> Lattice(dimension=2,randomize=False,order='natural').gen_samples(4, warn=False)
array([[0.  , 0.  ],
       [0.5 , 0.5 ],
       [0.25, 0.75],
       [0.75, 0.25]])
>>> Lattice(dimension=2,randomize=False,order='linear').gen_samples(4, warn=False)
array([[0.  , 0.  ],
       [0.25, 0.75],
       [0.5 , 0.5 ],
       [0.75, 0.25]])
>>> Lattice(dimension=2,randomize=False,order='mps').gen_samples(4, warn=False)
array([[0.  , 0.  ],
       [0.5 , 0.5 ],
       [0.25, 0.75],
       [0.75, 0.25]])

References

[1] Sou-Cheng T. Choi, Yuhan Ding, Fred J. Hickernell, Lan Jiang, Lluis Antoni Jimenez Rugama, Da Li, Jagadeeswaran Rathinavel, Xin Tong, Kan Zhang, Yizhi Zhang, and Xuan Zhou, GAIL: Guaranteed Automatic Integration Library (Version 2.3) [MATLAB Software], 2019. Available from http://gailgithub.github.io/GAIL_Dev/

[2] F.Y. Kuo & D. Nuyens. Application of quasi-Monte Carlo methods to elliptic PDEs with random diffusion coefficients - a survey of analysis and implementation, Foundations of Computational Mathematics, 16(6):1631-1696, 2016. springer link: https://link.springer.com/article/10.1007/s10208-016-9329-5 arxiv link: https://arxiv.org/abs/1606.06613

[3] D. Nuyens, The Magic Point Shop of QMC point generators and generating vectors. MATLAB and Python software, 2018. Available from https://people.cs.kuleuven.be/~dirk.nuyens/

[4] Constructing embedded lattice rules for multivariate integration R Cools, FY Kuo, D Nuyens - SIAM J. Sci. Comput., 28(6), 2162-2188.

[5] L’Ecuyer, Pierre & Munger, David. (2015). LatticeBuilder: A General Software Tool for Constructing Rank-1 Lattice Rules. ACM Transactions on Mathematical Software. 42. 10.1145/2754929.

__init__(dimension=1, randomize=True, order='natural', seed=None, generating_vector='lattice_vec.3600.20.npy', d_max=None, m_max=None)
Parameters

dimension (int or ndarray) –

dimension of the generator.

If an int is passed in, use sequence dimensions [0,…,dimensions-1]. If a ndarray is passed in, use these dimension indices in the sequence.

randomize (bool): If True, apply shift to generated samples. Note: Non-randomized lattice sequence includes the origin. order (str): ‘linear’, ‘natural’, or ‘mps’ ordering. seed (None or int or numpy.random.SeedSeq): seed the random number generator for reproducibility generating_vector (ndarray or str): generating matrix or path to generating matricies.

ndarray should have shape (d_max). a string generating_vector should be formatted like ‘lattice_vec.3600.20.npy’ where ‘name.d_max.m_max.npy’

d_max (int): maximum dimension m_max (int): 2^m_max is the max number of supported samples

Note

d_max and m_max are required if generating_vector is a ndarray. If generating_vector is an string (path), d_max and m_max can be taken from the file name if None

gen_samples(n=None, n_min=0, n_max=8, warn=True, return_unrandomized=False)

Generate lattice samples

Parameters
  • n (int) – if n is supplied, generate from n_min=0 to n_max=n samples. Otherwise use the n_min and n_max explicitly supplied as the following 2 arguments

  • n_min (int) – Starting index of sequence.

  • n_max (int) – Final index of sequence.

  • return_unrandomized (bool) – return samples without randomization as 2nd return value. Will not be returned if randomize=False.

Returns

(n_max-n_min) x d (dimension) array of samples

Return type

ndarray

Note

Lattice generates in blocks from 2**m to 2**(m+1) so generating n_min=3 to n_max=9 requires necessarily produces samples from n_min=2 to n_max=16 and automatically subsets. May be inefficient for non-powers-of-2 samples sizes.

pdf(x)

pdf of a standard uniform

Halton

class qmcpy.discrete_distribution.halton.Halton(dimension=1, randomize=True, generalize=True, seed=None)

Quasi-Random Halton nets.

>>> h_qrng = Halton(2,randomize='QRNG',generalize=True,seed=7)
>>> h_qrng.gen_samples(4)
array([[0.35362988, 0.38733489],
       [0.85362988, 0.72066823],
       [0.10362988, 0.05400156],
       [0.60362988, 0.498446  ]])
>>> h_qrng.gen_samples(1)
array([[0.35362988, 0.38733489]])
>>> h_qrng
Halton (DiscreteDistribution Object)
    d               2^(1)
    dvec            [0 1]
    randomize       QRNG
    generalize      1
    entropy         7
    spawn_key       ()
>>> h_owen = Halton(2,randomize='OWEN',generalize=False,seed=7)
>>> h_owen.gen_samples(4)
array([[0.64637012, 0.48226667],
       [0.14637012, 0.81560001],
       [0.89637012, 0.14893334],
       [0.39637012, 0.59337779]])
>>> h_owen
Halton (DiscreteDistribution Object)
    d               2^(1)
    dvec            [0 1]
    randomize       OWEN
    generalize      0
    entropy         7
    spawn_key       ()

References

[1] Marius Hofert and Christiane Lemieux (2019). qrng: (Randomized) Quasi-Random Number Generators. R package version 0.0-7. https://CRAN.R-project.org/package=qrng.

[2] Owen, A. B. “A randomized Halton algorithm in R,” 2017. arXiv:1706.02808 [stat.CO]

__init__(dimension=1, randomize=True, generalize=True, seed=None)
Parameters
  • dimension (int or ndarray) – dimension of the generator. If an int is passed in, use sequence dimensions [0,…,dimensions-1]. If a ndarray is passed in, use these dimension indices in the sequence.

  • randomize (str/bool) – select randomization method from ‘QRNG” [1], (max dimension = 360, supports generalize=True, default if randomize=True) or ‘OWEN’ [2], (max dimension = 1000),

  • generalize (bool) – generalize flag, only applicable to the QRNG generator

  • seed (None or int or numpy.random.SeedSeq) – seed the random number generator for reproducibility

gen_samples(n=None, n_min=0, n_max=8, warn=True)

Generate samples

Parameters
  • n (int) – if n is supplied, generate from n_min=0 to n_max=n samples. Otherwise use the n_min and n_max explicitly supplied as the following 2 arguments

  • n_min (int) – Starting index of sequence.

  • n_max (int) – Final index of sequence.

Returns

(n_max-n_min) x d (dimension) array of samples

Return type

ndarray

halton_owen(n, n0, d0=0)

see gen_samples method and [1] Owen, A. B.A randomized Halton algorithm in R2017. arXiv:1706.02808 [stat.CO].

pdf(x)

ABSTRACT METHOD to evaluate pdf of distribution the samples mimic at locations of x.

IID Standard Uniform

class qmcpy.discrete_distribution.iid_std_uniform.IIDStdUniform(dimension=1, seed=None)

A wrapper around NumPy’s IID Standard Uniform generator numpy.random.rand.

>>> dd = IIDStdUniform(dimension=2,seed=7)
>>> dd.gen_samples(4)
array([[0.04386058, 0.58727432],
       [0.3691824 , 0.65212985],
       [0.69669968, 0.10605352],
       [0.63025643, 0.13630282]])
>>> dd
IIDStdUniform (DiscreteDistribution Object)
    d               2^(1)
    entropy         7
    spawn_key       ()
__init__(dimension=1, seed=None)
Parameters
  • dimension (int) – dimension of samples

  • seed (None or int or numpy.random.SeedSeq) – seed the random number generator for reproducibility

gen_samples(n)

Generate samples

Parameters

n (int) – Number of observations to generate

Returns

n x self.d array of samples

Return type

ndarray

pdf(x)

ABSTRACT METHOD to evaluate pdf of distribution the samples mimic at locations of x.

True Measure Class

_images/true_measure_specific.png

Abstract Measure Class

class qmcpy.true_measure._true_measure.TrueMeasure

True Measure abstract class. DO NOT INSTANTIATE.

gen_samples(*args, **kwargs)

Generate samples from the discrete distribution and transform them via the transform method.

Parameters
  • args (tuple) – positional arguments to the discrete distributions gen_samples method

  • kwargs (dict) – keyword arguments to the discrete distributions gen_samples method

Returns

n x d matrix of transformed samples

Return type

ndarray

spawn(s=1, dimensions=None)

Spawn new instances of the current discrete distribution but with new seeds and dimensions. Developed for multi-level and multi-replication (Q)MC algorithms.

Parameters
  • s (int) – number of spawn

  • dimensions (ndarray) – length s array of dimension for each spawn. Defaults to current dimension

Returns

list of TrueMeasures linked to newly spawned DiscreteDistributions

Return type

list

Uniform

class qmcpy.true_measure.uniform.Uniform(sampler, lower_bound=0.0, upper_bound=1.0)
>>> u = Uniform(DigitalNetB2(2,seed=7),lower_bound=[0,.5],upper_bound=[2,3])
>>> u.gen_samples(4)
array([[1.12538017, 0.93444992],
       [0.693306  , 2.12676579],
       [1.64149095, 2.88726434],
       [0.20844522, 1.73645241]])
>>> u
Uniform (TrueMeasure Object)
    lower_bound     [0.  0.5]
    upper_bound     [2 3]
__init__(sampler, lower_bound=0.0, upper_bound=1.0)
Parameters
  • sampler (DiscreteDistribution/TrueMeasure) – A discrete distribution from which to transform samples or a true measure by which to compose a transform

  • lower_bound (float) – a for Uniform(a,b)

  • upper_bound (float) – b for Uniform(a,b)

Gaussian

class qmcpy.true_measure.gaussian.Gaussian(sampler, mean=0.0, covariance=1.0, decomp_type='PCA')

Normal Measure.

>>> g = Gaussian(DigitalNetB2(2,seed=7),mean=[1,2],covariance=[[9,4],[4,5]])
>>> g.gen_samples(4)
array([[-0.23979685,  2.98944192],
       [ 2.45994002,  2.17853622],
       [-0.22923897, -1.92667105],
       [ 4.6127697 ,  4.25820377]])
>>> g
Gaussian (TrueMeasure Object)
    mean            [1 2]
    covariance      [[9 4]
                    [4 5]]
    decomp_type     PCA
__init__(sampler, mean=0.0, covariance=1.0, decomp_type='PCA')
Parameters
  • sampler (DiscreteDistribution/TrueMeasure) – A discrete distribution from which to transform samples or a true measure by which to compose a transform

  • mean (float) – mu for Normal(mu,sigma^2)

  • covariance (ndarray) – sigma^2 for Normal(mu,sigma^2). A float or d (dimension) vector input will be extended to covariance*eye(d)

  • decomp_type (str) – method of decomposition either “PCA” for principal component analysis or “Cholesky” for cholesky decomposition.

class qmcpy.true_measure.gaussian.Normal(sampler, mean=0.0, covariance=1.0, decomp_type='PCA')

Brownian Motion

class qmcpy.true_measure.brownian_motion.BrownianMotion(sampler, t_final=1, drift=0, decomp_type='PCA')

Geometric Brownian Motion.

>>> bm = BrownianMotion(DigitalNetB2(4,seed=7),t_final=2,drift=2)
>>> bm.gen_samples(2)
array([[0.44018403, 1.62690477, 2.69418273, 4.21753174],
       [1.97549563, 2.27002956, 2.92802765, 4.77126959]])
>>> bm
BrownianMotion (TrueMeasure Object)
    time_vec        [0.5 1.  1.5 2. ]
    drift           2^(1)
    mean            [1. 2. 3. 4.]
    covariance      [[0.5 0.5 0.5 0.5]
                    [0.5 1.  1.  1. ]
                    [0.5 1.  1.5 1.5]
                    [0.5 1.  1.5 2. ]]
    decomp_type     PCA
__init__(sampler, t_final=1, drift=0, decomp_type='PCA')
Parameters
  • sampler (DiscreteDistribution/TrueMeasure) – A discrete distribution from which to transform samples or a true measure by which to compose a transform

  • t_final (float) – end time for the Brownian Motion.

  • drift (int) – Gaussian mean is time_vec*drift

  • decomp_type (str) – method of decomposition either “PCA” for principal component analysis or “Cholesky” for cholesky decomposition.

Lebesgue

class qmcpy.true_measure.lebesgue.Lebesgue(sampler)
>>> Lebesgue(Gaussian(DigitalNetB2(2,seed=7)))
Lebesgue (TrueMeasure Object)
    transform       Gaussian (TrueMeasure Object)
                       mean            0
                       covariance      1
                       decomp_type     PCA
>>> Lebesgue(Uniform(DigitalNetB2(2,seed=7)))
Lebesgue (TrueMeasure Object)
    transform       Uniform (TrueMeasure Object)
                       lower_bound     0
                       upper_bound     1
__init__(sampler)
Parameters

sampler (TrueMeasure) – A true measure by which to compose a transform.

Continuous Bernoulli

class qmcpy.true_measure.bernoulli_cont.BernoulliCont(sampler, lam=0.5)
>>> bc = BernoulliCont(DigitalNetB2(2,seed=7),lam=.2)
>>> bc.gen_samples(4)
array([[0.39545122, 0.10073414],
       [0.21719142, 0.48293404],
       [0.68958314, 0.90847415],
       [0.05871131, 0.33436033]])
>>> bc
BernoulliCont (TrueMeasure Object)
    lam             0.200

See https://en.wikipedia.org/wiki/Continuous_Bernoulli_distribution

__init__(sampler, lam=0.5)
Parameters
  • sampler (DiscreteDistribution/TrueMeasure) – A discrete distribution from which to transform samples or a true measure by which to compose a transform

  • lam (ndarray) – 0 < lambda < 1, a shape parameter, independent for each dimension

Johnson’s SU

class qmcpy.true_measure.johnsons_su.JohnsonsSU(sampler, gamma=1, xi=1, delta=2, lam=2)
>>> jsu = JohnsonsSU(DigitalNetB2(2,seed=7),gamma=1,xi=2,delta=3,lam=4)
>>> jsu.gen_samples(4)
array([[ 0.86224892, -0.76967276],
       [ 0.07317047,  1.17727769],
       [ 1.89093286,  2.9341619 ],
       [-1.30283298,  0.62269632]])
>>> jsu
JohnsonsSU (TrueMeasure Object)
    gamma           1
    xi              2^(1)
    delta           3
    lam             2^(2)

See https://en.wikipedia.org/wiki/Johnson%27s_SU-distribution

__init__(sampler, gamma=1, xi=1, delta=2, lam=2)
Parameters
  • sampler (DiscreteDistribution/TrueMeasure) – A discrete distribution from which to transform samples or a true measure by which to compose a transform

  • gamma (ndarray) – gamma

  • xi (ndarray) – xi

  • delta (ndarray) – delta > 0

  • lam (ndarray) – lambda > 0

Kumaraswamy

class qmcpy.true_measure.kumaraswamy.Kumaraswamy(sampler, a=2, b=2)
>>> k = Kumaraswamy(DigitalNetB2(2,seed=7),a=[1,2],b=[3,4])
>>> k.gen_samples(4)
array([[0.24096272, 0.21587652],
       [0.13227662, 0.4808615 ],
       [0.43615893, 0.73428949],
       [0.03602294, 0.39602319]])
>>> k
Kumaraswamy (TrueMeasure Object)
    a               [1 2]
    b               [3 4]

See https://en.wikipedia.org/wiki/Kumaraswamy_distribution

__init__(sampler, a=2, b=2)
Parameters
  • sampler (DiscreteDistribution/TrueMeasure) – A discrete distribution from which to transform samples or a true measure by which to compose a transform

  • a (ndarray) – alpha > 0

  • b (ndarray) – beta > 0

SciPy Wrapper

class qmcpy.true_measure.scipy_wrapper.SciPyWrapper(sampler, scipy_distrib, **scipy_distrib_kwargs)
>>> triangular = SciPyWrapper(
...     sampler = DigitalNetB2(2,seed=7),
...     scipy_distrib = scipy.stats.triang,
...     c = [0.1,.2],
...     loc = [1,2],
...     scale = [3,4])
>>> triangular.gen_samples(4)
array([[2.11792393, 2.74571838],
       [1.6995412 , 3.88553573],
       [2.79502629, 5.24025887],
       [1.30634136, 3.45650562]])
>>> triangular
SciPyWrapper (TrueMeasure Object)
    scipy_distrib   triang
    c               [0.1 0.2]
    loc             [1 2]
    scale           [3 4]
__init__(sampler, scipy_distrib, **scipy_distrib_kwargs)
Parameters

Integrand Class

_images/integrand_specific.png

Abstract Integrand Class

class qmcpy.integrand._integrand.Integrand

Integrand abstract class. DO NOT INSTANTIATE.

bound_fun(bound_low, bound_high)

compute the bounds on the combined function based on bounds for the individual functions. Defaults to the identity where we essentiallly do not combine integrands, but instead integrate each function individually.

Parameters
  • bound_low (ndarray) – length Integrand.dprime lower error bound

  • bound_high (ndarray) – length Integrand.dprime upper error bound

Returns

lower bound on function combining estimates, ndarray: upper bound on function combining estimates, bool ndarray: flags to override suffient combined integrand estimation.

e.g. when approximating a ratio of integrals, if the denominator’s bounds straddle 0, then returning True here forces ratio to be flagged as insuffiently approximated.

Return type

ndarray

dependency(flags_comb)

takes a vector of indicators of wheather of not the error bound is satisfied for combined integrands and which returns flags for individual integrands. For example, if we are taking the ratio of 2 individual integrands, then getting flag_comb=True means the ratio has not been approximated to within the tolerance, so the dependency function should return [True,True] indicating that both the numerator and denominator integrands need to be better approximated. :param flags_comb: flags indicating wheather the combined integrals are insufficiently approximated :type flags_comb: bool ndarray

Returns

length (Integrand.dprime) flags for individual integrands

Return type

(bool ndarray)

f(x, periodization_transform='NONE', *args, **kwargs)

Evalute transformed integrand based on true measures and discrete distribution

Parameters
  • x (ndarray) – n x d array of samples from a discrete distribution

  • ptransform (str) – periodization transform.

  • *args – other ordered args to g

  • **kwargs (dict) – other keyword args to g

Returns

length n vector of funciton evaluations

Return type

ndarray

g(t, *args, **kwargs)

ABSTRACT METHOD for original integrand to be integrated.

Parameters

t (ndarray) – n x d array of samples to be intput into orignal integrand.

Returns

n vector of function evaluations

Return type

ndarray

spawn(levels)

Spawn new instances of the current integrand at the specified levels.

Parameters

levels (ndarray) – array of levels at which to spawn new integrands

Returns

list of Integrands linked to newly spawned TrueMeasures and DiscreteDistributions

Return type

list

Custom Function

class qmcpy.integrand.custom_fun.CustomFun(true_measure, g, dprime=1)

Integrand wrapper for a user’s function

>>> cf = CustomFun(
...     true_measure = Gaussian(DigitalNetB2(2,seed=7),mean=[1,2]),
...     g = lambda x: x[:,0]**2*x[:,1])
>>> x = cf.discrete_distrib.gen_samples(2**10)
>>> y = cf.f(x)
>>> y.shape
(1024, 1)
>>> y.mean()
3.995...
>>> cf = CustomFun(
...     true_measure = Uniform(DigitalNetB2(3,seed=7),lower_bound=[2,3,4],upper_bound=[4,5,6]),
...     g = lambda x: x,
...     dprime = 3)
>>> x = cf.discrete_distrib.gen_samples(2**10)
>>> y = cf.f(x)
>>> y.shape
(1024, 3)
>>> y.mean(0)
array([3., 4., 5.])
__init__(true_measure, g, dprime=1)
Parameters
  • true_measure (TrueMeasure) – a TrueMeasure instance.

  • g (function) – a function handle.

  • dprime (int) – output dimension of the function.

g(t, *args, **kwargs)

ABSTRACT METHOD for original integrand to be integrated.

Parameters

t (ndarray) – n x d array of samples to be intput into orignal integrand.

Returns

n vector of function evaluations

Return type

ndarray

Keister Function

class qmcpy.integrand.keister.Keister(sampler)

f(\boldsymbol{t}) = \pi^{d/2} \cos(\| \boldsymbol{t} \|).

The standard example integrates the Keister integrand with respect to an IID Gaussian distribution with variance 1./2.

>>> k = Keister(DigitalNetB2(2,seed=7))
>>> x = k.discrete_distrib.gen_samples(2**10)
>>> y = k.f(x)
>>> y.mean()
1.808...
>>> k.true_measure
Gaussian (TrueMeasure Object)
    mean            0
    covariance      2^(-1)
    decomp_type     PCA
>>> k = Keister(Gaussian(DigitalNetB2(2,seed=7),mean=0,covariance=2))
>>> x = k.discrete_distrib.gen_samples(2**12)
>>> y = k.f(x)
>>> y.mean()
1.808...
>>> yp = k.f(x,periodization_transform='c2sin')
>>> yp.mean()
1.807...

References

[1] B. D. Keister, Multidimensional Quadrature Algorithms, Computers in Physics, 10, pp. 119-122, 1996.

__init__(sampler)
Parameters

sampler (DiscreteDistribution/TrueMeasure) – A discrete distribution from which to transform samples or a true measure by which to compose a transform

g(t)

ABSTRACT METHOD for original integrand to be integrated.

Parameters

t (ndarray) – n x d array of samples to be intput into orignal integrand.

Returns

n vector of function evaluations

Return type

ndarray

Box Integral

class qmcpy.integrand.box_integral.BoxIntegral(sampler, s=array([1, 2]))

B_s(x) = \left(\sum_{j=1}^d x_j^2 \right)^{s/2}

>>> l1 = BoxIntegral(DigitalNetB2(2,seed=7), s=[7])
>>> x1 = l1.discrete_distrib.gen_samples(2**10)
>>> y1 = l1.f(x1)
>>> y1.shape
(1024, 1)
>>> y1.mean(0)
array([0.75156724])
>>> l2 = BoxIntegral(DigitalNetB2(5,seed=7), s=[-7,7])
>>> x2 = l2.discrete_distrib.gen_samples(2**10)
>>> y2 = l2.f(x2,compute_flags=[1,1])
>>> y2.shape
(1024, 2)
>>> y2.mean(0)
array([ 6.67548708, 10.52267786])

References:

[1] D.H. Bailey, J.M. Borwein, R.E. Crandall,Box integrals, Journal of Computational and Applied Mathematics, Volume 206, Issue 1, 2007, Pages 196-208, ISSN 0377-0427, https://doi.org/10.1016/j.cam.2006.06.010. (https://www.sciencedirect.com/science/article/pii/S0377042706004250)

[2] https://www.davidhbailey.com/dhbpapers/boxintegrals.pdf

__init__(sampler, s=array([1, 2]))
Parameters
  • sampler (DiscreteDistribution/TrueMeasure) – A discrete distribution from which to transform samples or a true measure by which to compose a transform

  • s (list or ndarray) – vectorized s parameter, len(s) is the number of vectorized integrals to evalute.

g(t, **kwargs)

ABSTRACT METHOD for original integrand to be integrated.

Parameters

t (ndarray) – n x d array of samples to be intput into orignal integrand.

Returns

n vector of function evaluations

Return type

ndarray

European Option

class qmcpy.integrand.european_option.EuropeanOption(sampler, volatility=0.5, start_price=30, strike_price=35, interest_rate=0, t_final=1, call_put='call')

European financial option.

>>> eo = EuropeanOption(DigitalNetB2(4,seed=7),call_put='put')
>>> eo
EuropeanOption (Integrand Object)
    volatility      2^(-1)
    call_put        put
    start_price     30
    strike_price    35
    interest_rate   0
>>> x = eo.discrete_distrib.gen_samples(2**12)
>>> y = eo.f(x)
>>> y.mean()
9.209...
>>> eo = EuropeanOption(BrownianMotion(DigitalNetB2(4,seed=7),drift=1),call_put='put')
>>> x = eo.discrete_distrib.gen_samples(2**12)
>>> y = eo.f(x)
>>> y.mean()
9.162...
>>> eo.get_exact_value()
9.211452976234058
__init__(sampler, volatility=0.5, start_price=30, strike_price=35, interest_rate=0, t_final=1, call_put='call')
Parameters
  • sampler (DiscreteDistribution/TrueMeasure) – A discrete distribution from which to transform samples or a true measure by which to compose a transform

  • volatility (float) – sigma, the volatility of the asset

  • start_price (float) – S(0), the asset value at t=0

  • strike_price (float) – strike_price, the call/put offer

  • interest_rate (float) – r, the annual interest rate

  • t_final (float) – exercise time

  • call_put (str) – ‘call’ or ‘put’ option

g(t)

See abstract method.

get_exact_value()

Get the fair price of a European call/put option.

Returns

fair price

Return type

float

Asian Option

class qmcpy.integrand.asian_option.AsianOption(sampler, volatility=0.5, start_price=30.0, strike_price=35.0, interest_rate=0.0, t_final=1, call_put='call', mean_type='arithmetic', multilevel_dims=None, _dim_frac=0)

Asian financial option.

>>> ac = AsianOption(DigitalNetB2(4,seed=7))
>>> ac
AsianOption (Integrand Object)
    volatility      2^(-1)
    call_put        call
    start_price     30
    strike_price    35
    interest_rate   0
    mean_type       arithmetic
    dim_frac        0
>>> x = ac.discrete_distrib.gen_samples(2**12)
>>> y = ac.f(x)
>>> y.mean()
1.768...
>>> level_dims = [2,4,8]
>>> ac2_multilevel = AsianOption(DigitalNetB2(seed=7),multilevel_dims=level_dims)
>>> levels_to_spawn = arange(ac2_multilevel.max_level+1)
>>> ac2_single_levels = ac2_multilevel.spawn(levels_to_spawn)
>>> yml = 0
>>> for ac2_single_level in ac2_single_levels:
...     x = ac2_single_level.discrete_distrib.gen_samples(2**12)
...     level_est = ac2_single_level.f(x).mean()
...     yml += level_est
>>> yml
1.779...
__init__(sampler, volatility=0.5, start_price=30.0, strike_price=35.0, interest_rate=0.0, t_final=1, call_put='call', mean_type='arithmetic', multilevel_dims=None, _dim_frac=0)
Parameters
  • sampler (DiscreteDistribution/TrueMeasure) – A discrete distribution from which to transform samples or a true measure by which to compose a transform

  • volatility (float) – sigma, the volatility of the asset

  • start_price (float) – S(0), the asset value at t=0

  • strike_price (float) – strike_price, the call/put offer

  • interest_rate (float) – r, the annual interest rate

  • t_final (float) – exercise time

  • mean_type (string) – ‘arithmetic’ or ‘geometric’ mean

  • multilevel_dims (list of ints) – list of dimensions at each level. Leave as None for single-level problems

  • _dim_frac (float) – for internal use only, users should not set this parameter.

g(t)

ABSTRACT METHOD for original integrand to be integrated.

Parameters

t (ndarray) – n x d array of samples to be intput into orignal integrand.

Returns

n vector of function evaluations

Return type

ndarray

Multilevel Call Options with Milstein Discretization

class qmcpy.integrand.ml_call_options.MLCallOptions(sampler, option='european', volatility=0.2, start_strike_price=100.0, interest_rate=0.05, t_final=1.0, _level=0)

Various call options from finance using Milstein discretization with 2^l timesteps on level l.

>>> mlco_original = MLCallOptions(DigitalNetB2(seed=7))
>>> mlco_original
MLCallOptions (Integrand Object)
    option          european
    sigma           0.200
    k               100
    r               0.050
    t               1
    b               85
    level           0
>>> mlco_ml_dims = mlco_original.spawn(levels=arange(4))
>>> yml = 0
>>> for mlco in mlco_ml_dims:
...     x = mlco.discrete_distrib.gen_samples(2**10)
...     yml += mlco.f(x).mean()
>>> yml
10.393...

References:

[1] M.B. Giles. Improved multilevel Monte Carlo convergence using the Milstein scheme. 343-358, in Monte Carlo and Quasi-Monte Carlo Methods 2006, Springer, 2008. http://people.maths.ox.ac.uk/~gilesm/files/mcqmc06.pdf.

__init__(sampler, option='european', volatility=0.2, start_strike_price=100.0, interest_rate=0.05, t_final=1.0, _level=0)
Parameters
  • sampler (DiscreteDistribution/TrueMeasure) – A discrete distribution from which to transform samples or a true measure by which to compose a transform

  • option_type (str) – type of option in [“European”,”Asian”]

  • volatility (float) – sigma, the volatility of the asset

  • start_strike_price (float) – S(0), the asset value at t=0, and K, the strike price. Assume start_price = strike_price

  • interest_rate (float) – r, the annual interest rate

  • t_final (float) – exercise time

  • _level (int) – for internal use only, users should not set this parameter.

g(t)
Parameters

t (ndarray) – Gaussian(0,1^2) samples

Returns

First, an ndarray of length 6 vector of summary statistic sums. Second, a float of cost on this level.

Return type

tuple

get_exact_value()

Print exact analytic value, based on s0=k.

Linear Function

class qmcpy.integrand.linear0.Linear0(sampler)
>>> l = Linear0(DigitalNetB2(100,seed=7))
>>> x = l.discrete_distrib.gen_samples(2**10)
>>> y = l.f(x)
>>> y.mean()
-1.175...e-08
>>> ytf = l.f(x,periodization_transform='C1SIN')
>>> ytf.mean()
-4.050...e-12
__init__(sampler)
Parameters

sampler (DiscreteDistribution/TrueMeasure) – A discrete distribution from which to transform samples or a true measure by which to compose a transform

g(t)

ABSTRACT METHOD for original integrand to be integrated.

Parameters

t (ndarray) – n x d array of samples to be intput into orignal integrand.

Returns

n vector of function evaluations

Return type

ndarray

Sobol’ Indices

class qmcpy.integrand.sobol_indices.SobolIndices(integrand, indices='singletons')

Sobol’ Indicies in QMCPy.

>>> dnb2 = DigitalNetB2(dimension=3,seed=7)
>>> keister_d = Keister(dnb2)
>>> keister_indices = SobolIndices(keister_d,indices='singletons')
>>> sc = CubQMCCLT(keister_indices,abs_tol=.05)
>>> solution,data = sc.integrate()
>>> data
MeanVarDataRep (AccumulateData Object)
    solution        [0.321 0.331 0.328 0.334 0.341 0.334]
    indv_error_bound [0.092 0.103 0.046 0.048 0.072 0.058 0.007 0.051]
    ci_low          [1.555 1.595 1.64  1.667 1.679 1.659 2.163 9.77 ]
    ci_high         [1.739 1.801 1.731 1.763 1.823 1.775 2.177 9.871]
    ci_comb_low     [0.3   0.307 0.316 0.321 0.324 0.32 ]
    ci_comb_high    [0.342 0.354 0.34  0.346 0.358 0.349]
    solution_comb   [0.321 0.331 0.328 0.334 0.341 0.334]
    flags_comb      [False False False False False False]
    flags_indv      [False False False False False False False False]
    n_total         2^(12)
    n               [256. 256. 256. 256. 256. 256. 256. 256.]
    replications    2^(4)
    time_integrate  ...
CubQMCCLT (StoppingCriterion Object)
    inflate         1.200
    alpha           0.010
    abs_tol         0.050
    rel_tol         0
    n_init          2^(8)
    n_max           2^(30)
SobolIndices (Integrand Object)
    indices         [[0]
                    [1]
                    [2]]
Gaussian (TrueMeasure Object)
    mean            0
    covariance      2^(-1)
    decomp_type     PCA
DigitalNetB2 (DiscreteDistribution Object)
    d               6
    dvec            [0 1 2 3 4 5]
    randomize       LMS_DS
    graycode        0
    entropy         7
    spawn_key       (0,)

References

[1] Art B. Owen.Monte Carlo theory, methods and examples. 2013. Appendix A.

__init__(integrand, indices='singletons')
Parameters
  • integrand (Integrand) – integrand to find Sobol’ indices of

  • indices (list of lists) – each element of indices should be a list of indices, u, at which to compute the Sobol’ indices. The default indices=’singletons’ sets indices=[[0],[1],…[d-1]]. Should not include [], the null set

bound_fun(bound_low, bound_high)

compute the bounds on the combined function based on bounds for the individual functions. Defaults to the identity where we essentiallly do not combine integrands, but instead integrate each function individually.

Parameters
  • bound_low (ndarray) – length Integrand.dprime lower error bound

  • bound_high (ndarray) – length Integrand.dprime upper error bound

Returns

lower bound on function combining estimates, ndarray: upper bound on function combining estimates, bool ndarray: flags to override suffient combined integrand estimation.

e.g. when approximating a ratio of integrals, if the denominator’s bounds straddle 0, then returning True here forces ratio to be flagged as insuffiently approximated.

Return type

ndarray

dependency(flags_comb)

takes a vector of indicators of wheather of not the error bound is satisfied for combined integrands and which returns flags for individual integrands. For example, if we are taking the ratio of 2 individual integrands, then getting flag_comb=True means the ratio has not been approximated to within the tolerance, so the dependency function should return [True,True] indicating that both the numerator and denominator integrands need to be better approximated. :param flags_comb: flags indicating wheather the combined integrals are insufficiently approximated :type flags_comb: bool ndarray

Returns

length (Integrand.dprime) flags for individual integrands

Return type

(bool ndarray)

f(x, periodization_transform='NONE', *args, **kwargs)

Evalute transformed integrand based on true measures and discrete distribution

Parameters
  • x (ndarray) – n x d array of samples from a discrete distribution

  • ptransform (str) – periodization transform.

  • *args – other ordered args to g

  • **kwargs (dict) – other keyword args to g

Returns

length n vector of funciton evaluations

Return type

ndarray

Stopping Criterion Algorithms

_images/stopping_criterion_specific.png

Abstract Stopping Criterion Class

class qmcpy.stopping_criterion._stopping_criterion.StoppingCriterion(allowed_levels, allowed_distribs, allow_vectorized_integrals)

Stopping Criterion abstract class. DO NOT INSTANTIATE.

__init__(allowed_levels, allowed_distribs, allow_vectorized_integrals)
Parameters
  • distribution (DiscreteDistribution) – a DiscreteDistribution

  • allowed_levels (list) – which integrand types are supported: ‘single’, ‘fixed-multi’, ‘adaptive-multi’

  • allowed_distribs (list) – list of compatible DiscreteDistribution classes

integrate()

ABSTRACT METHOD to determine the number of samples needed to satisfy the tolerance.

Returns

tuple containing:
  • solution (float): approximation to the integral

  • data (AccumulateData): an AccumulateData object

Return type

tuple

set_tolerance(*args, **kwargs)

ABSTRACT METHOD to reset the absolute tolerance.

Guaranteed Digital Net Cubature (QMC)

class qmcpy.stopping_criterion.cub_qmc_net_g.CubQMCNetG(integrand, abs_tol=0.01, rel_tol=0.0, n_init=1024.0, n_max=34359738368.0, fudge=<function CubQMCNetG.<lambda>>, check_cone=False, control_variates=[], control_variate_means=[], update_beta=False)

Quasi-Monte Carlo method using Sobol’ cubature over the d-dimensional region to integrate within a specified generalized error tolerance with guarantees under Walsh-Fourier coefficients cone decay assumptions.

>>> k = Keister(DigitalNetB2(2,seed=7))
>>> sc = CubQMCNetG(k,abs_tol=.05)
>>> solution,data = sc.integrate()
>>> data
LDTransformData (AccumulateData Object)
    solution        1.809
    error_bound     0.005
    n_total         2^(10)
    time_integrate  ...
CubQMCNetG (StoppingCriterion Object)
    abs_tol         0.050
    rel_tol         0
    n_init          2^(10)
    n_max           2^(35)
Keister (Integrand Object)
Gaussian (TrueMeasure Object)
    mean            0
    covariance      2^(-1)
    decomp_type     PCA
DigitalNetB2 (DiscreteDistribution Object)
    d               2^(1)
    dvec            [0 1]
    randomize       LMS_DS
    graycode        0
    entropy         7
    spawn_key       ()
>>> dd = DigitalNetB2(3,seed=7)
>>> g1 = CustomFun(Uniform(dd,0,2),lambda t: 10*t[:,0]-5*t[:,1]**2+t[:,2]**3)
>>> cv1 = CustomFun(Uniform(dd,0,2),lambda t: t[:,0])
>>> cv2 = CustomFun(Uniform(dd,0,2),lambda t: t[:,1]**2)
>>> sc = CubQMCNetG(g1,abs_tol=1e-6,check_cone=True,
...     control_variates = [cv1,cv2],
...     control_variate_means = [1,4/3])
>>> sol,data = sc.integrate()
>>> print(sol)
5.333...
>>> exactsol = 16/3
>>> print(abs(sol-exactsol)<1e-6)
True

Original Implementation:

References

[1] Fred J. Hickernell and Lluis Antoni Jimenez Rugama, Reliable adaptive cubature using digital sequences, 2014. Submitted for publication: arXiv:1410.8615.

[2] Sou-Cheng T. Choi, Yuhan Ding, Fred J. Hickernell, Lan Jiang, Lluis Antoni Jimenez Rugama, Da Li, Jagadeeswaran Rathinavel, Xin Tong, Kan Zhang, Yizhi Zhang, and Xuan Zhou, GAIL: Guaranteed Automatic Integration Library (Version 2.3) [MATLAB Software], 2019. Available from http://gailgithub.github.io/GAIL_Dev/

Guarantee:

This algorithm computes the integral of real valued functions in [0,1]^d with a prescribed generalized error tolerance. The Fourier coefficients of the integrand are assumed to be absolutely convergent. If the algorithm terminates without warning messages, the output is given with guarantees under the assumption that the integrand lies inside a cone of functions. The guarantee is based on the decay rate of the Fourier coefficients. For integration over domains other than [0,1]^d, this cone condition applies to f \circ \psi (the composition of the functions) where \psi is the transformation function for [0,1]^d to the desired region. For more details on how the cone is defined, please refer to the references below.

__init__(integrand, abs_tol=0.01, rel_tol=0.0, n_init=1024.0, n_max=34359738368.0, fudge=<function CubQMCNetG.<lambda>>, check_cone=False, control_variates=[], control_variate_means=[], update_beta=False)
Parameters
  • integrand (Integrand) – an instance of Integrand

  • abs_tol (float) – absolute error tolerance

  • rel_tol (float) – relative error tolerance

  • n_init (int) – initial number of samples

  • n_max (int) – maximum number of samples

  • fudge (function) – positive function multiplying the finite sum of Fast Fourier coefficients specified in the cone of functions

  • check_cone (boolean) – check if the function falls in the cone

  • control_variates (list) – list of integrand objects to be used as control variates. Control variates are currently only compatible with single level problems. The same discrete distribution instance must be used for the integrand and each of the control variates.

  • control_variate_means (list) – list of means for each control variate

  • update_beta (bool) – update control variate beta coefficients at each iteration?

class qmcpy.stopping_criterion.cub_qmc_net_g.CubQMCSobolG(integrand, abs_tol=0.01, rel_tol=0.0, n_init=1024.0, n_max=34359738368.0, fudge=<function CubQMCNetG.<lambda>>, check_cone=False, control_variates=[], control_variate_means=[], update_beta=False)

Guaranteed Lattice Cubature (QMC)

class qmcpy.stopping_criterion.cub_qmc_lattice_g.CubQMCLatticeG(integrand, abs_tol=0.01, rel_tol=0.0, n_init=1024.0, n_max=34359738368.0, fudge=<function CubQMCLatticeG.<lambda>>, check_cone=False, ptransform='Baker')

Stopping Criterion quasi-Monte Carlo method using rank-1 Lattices cubature over a d-dimensional region to integrate within a specified generalized error tolerance with guarantees under Fourier coefficients cone decay assumptions.

>>> k = Keister(Lattice(2,seed=7))
>>> sc = CubQMCLatticeG(k,abs_tol=.05)
>>> solution,data = sc.integrate()
>>> data
LDTransformData (AccumulateData Object)
    solution        1.810
    error_bound     0.005
    n_total         2^(10)
    time_integrate  ...
CubQMCLatticeG (StoppingCriterion Object)
    abs_tol         0.050
    rel_tol         0
    n_init          2^(10)
    n_max           2^(35)
Keister (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         7
    spawn_key       ()

Original Implementation:

References

[1] Lluis Antoni Jimenez Rugama and Fred J. Hickernell, “Adaptive multidimensional integration based on rank-1 lattices,” Monte Carlo and Quasi-Monte Carlo Methods: MCQMC, Leuven, Belgium, April 2014 (R. Cools and D. Nuyens, eds.), Springer Proceedings in Mathematics and Statistics, vol. 163, Springer-Verlag, Berlin, 2016, arXiv:1411.1966, pp. 407-422.

[2] Sou-Cheng T. Choi, Yuhan Ding, Fred J. Hickernell, Lan Jiang, Lluis Antoni Jimenez Rugama, Da Li, Jagadeeswaran Rathinavel, Xin Tong, Kan Zhang, Yizhi Zhang, and Xuan Zhou, GAIL: Guaranteed Automatic Integration Library (Version 2.3) [MATLAB Software], 2019. Available from http://gailgithub.github.io/GAIL_Dev/

Guarantee

This algorithm computes the integral of real valued functions in [0,1]^d with a prescribed generalized error tolerance. The Fourier coefficients of the integrand are assumed to be absolutely convergent. If the algorithm terminates without warning messages, the output is given with guarantees under the assumption that the integrand lies inside a cone of functions. The guarantee is based on the decay rate of the Fourier coefficients. For integration over domains other than [0,1]^d, this cone condition applies to f \circ \psi (the composition of the functions) where \psi is the transformation function for [0,1]^d to the desired region. For more details on how the cone is defined, please refer to the references below.

__init__(integrand, abs_tol=0.01, rel_tol=0.0, n_init=1024.0, n_max=34359738368.0, fudge=<function CubQMCLatticeG.<lambda>>, check_cone=False, ptransform='Baker')
Parameters
  • integrand (Integrand) – an instance of Integrand

  • abs_tol (float) – absolute error tolerance

  • rel_tol (float) – relative error tolerance

  • n_init (int) – initial number of samples

  • n_max (int) – maximum number of samples

  • fudge (function) – positive function multiplying the finite sum of Fast Fourier coefficients specified in the cone of functions

  • check_cone (boolean) – check if the function falls in the cone

Bayesian Lattice Cubature (QMC)

class qmcpy.stopping_criterion.cub_qmc_bayes_lattice_g.CubBayesLatticeG(integrand, abs_tol=0.01, rel_tol=0, n_init=256, n_max=4194304, order=2, alpha=0.01, ptransform='C1sin')

Stopping criterion for Bayesian Cubature using rank-1 Lattice sequence with guaranteed accuracy over a d-dimensional region to integrate within a specified generalized error tolerance with guarantees under Bayesian assumptions.

>>> k = Keister(Lattice(2, order='linear', seed=123456789))
>>> sc = CubBayesLatticeG(k,abs_tol=.05)
>>> solution,data = sc.integrate()
>>> data
LDTransformBayesData (AccumulateData Object)
    solution        1.808
    error_bound     7.37e-04
    n_total         256
    time_integrate  ...
CubBayesLatticeG (StoppingCriterion Object)
    abs_tol         0.050
    rel_tol         0
    n_init          2^(8)
    n_max           2^(22)
    order           2^(1)
Keister (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           linear
    entropy         123456789
    spawn_key       ()
Adapted from

https://github.com/GailGithub/GAIL_Dev/blob/master/Algorithms/IntegrationExpectation/cubBayesLattice_g.m

Reference

[1] Sou-Cheng T. Choi, Yuhan Ding, Fred J. Hickernell, Lan Jiang, Lluis Antoni Jimenez Rugama, Da Li, Jagadeeswaran Rathinavel, Xin Tong, Kan Zhang, Yizhi Zhang, and Xuan Zhou, GAIL: Guaranteed Automatic Integration Library (Version 2.3) [MATLAB Software], 2019. Available from http://gailgithub.github.io/GAIL_Dev/

Guarantee

This algorithm attempts to calculate the integral of function f over the hyperbox [0,1]^d to a prescribed error tolerance tolfun:= max(abstol,reltol*| I |) with guaranteed confidence level, e.g., 99% when alpha=0.5%. If the algorithm terminates without showing any warning messages and provides an answer Q, then the following inequality would be satisfied:

Pr(| Q - I | <= tolfun) = 99%

This Bayesian cubature algorithm guarantees for integrands that are considered to be an instance of a gaussian process that fall in the middle of samples space spanned. Where The sample space is spanned by the covariance kernel parametrized by the scale and shape parameter inferred from the sampled values of the integrand. For more details on how the covariance kernels are defined and the parameters are obtained, please refer to the references below.

integrate()

ABSTRACT METHOD to determine the number of samples needed to satisfy the tolerance.

Returns

tuple containing:
  • solution (float): approximation to the integral

  • data (AccumulateData): an AccumulateData object

Return type

tuple

Bayesian Digital Net Cubature (QMC)

class qmcpy.stopping_criterion.cub_qmc_bayes_net_g.CubBayesNetG(integrand, abs_tol=0.01, rel_tol=0, n_init=256, n_max=4194304, alpha=0.01)

Stopping criterion for Bayesian Cubature using digital net sequence with guaranteed accuracy over a d-dimensional region to integrate within a specified generalized error tolerance with guarantees under Bayesian assumptions.

>>> k = Keister(DigitalNetB2(2, seed=123456789))
>>> sc = CubBayesNetG(k,abs_tol=.05)
>>> solution,data = sc.integrate()
>>> data
LDTransformBayesData (AccumulateData Object)
    solution        1.812
    error_bound     0.015
    n_total         256
    time_integrate  ...
CubBayesNetG (StoppingCriterion Object)
    abs_tol         0.050
    rel_tol         0
    n_init          2^(8)
    n_max           2^(22)
Keister (Integrand Object)
Gaussian (TrueMeasure Object)
    mean            0
    covariance      2^(-1)
    decomp_type     PCA
DigitalNetB2 (DiscreteDistribution Object)
    d               2^(1)
    dvec            [0 1]
    randomize       LMS_DS
    graycode        0
    entropy         123456789
    spawn_key       ()
Adapted from

https://github.com/GailGithub/GAIL_Dev/blob/master/Algorithms/IntegrationExpectation/cubBayesNet_g.m

Reference

[1] Sou-Cheng T. Choi, Yuhan Ding, Fred J. Hickernell, Lan Jiang, Lluis Antoni Jimenez Rugama, Da Li, Jagadeeswaran Rathinavel, Xin Tong, Kan Zhang, Yizhi Zhang, and Xuan Zhou, GAIL: Guaranteed Automatic Integration Library (Version 2.3) [MATLAB Software], 2019. Available from http://gailgithub.github.io/GAIL_Dev/

Guarantee

This algorithm attempts to calculate the integral of function f over the hyperbox [0,1]^d to a prescribed error tolerance tolfun:= max(abstol,reltol*| I |) with guaranteed confidence level, e.g., 99% when alpha=0.5%. If the algorithm terminates without showing any warning messages and provides an answer Q, then the following inequality would be satisfied:

Pr(| Q - I | <= tolfun) = 99%

This Bayesian cubature algorithm guarantees for integrands that are considered to be an instance of a gaussian process that fall in the middle of samples space spanned. Where The sample space is spanned by the covariance kernel parametrized by the scale and shape parameter inferred from the sampled values of the integrand. For more details on how the covariance kernels are defined and the parameters are obtained, please refer to the references below.

integrate()

ABSTRACT METHOD to determine the number of samples needed to satisfy the tolerance.

Returns

tuple containing:
  • solution (float): approximation to the integral

  • data (AccumulateData): an AccumulateData object

Return type

tuple

class qmcpy.stopping_criterion.cub_qmc_bayes_net_g.CubBayesSobolG(integrand, abs_tol=0.01, rel_tol=0, n_init=256, n_max=4194304, alpha=0.01)

CLT QMC Cubature (with Replications)

class qmcpy.stopping_criterion.cub_qmc_clt.CubQMCCLT(integrand, abs_tol=0.01, rel_tol=0.0, n_init=256.0, n_max=1073741824, inflate=1.2, alpha=0.01, replications=16.0, error_fun=<function CubQMCCLT.<lambda>>)

Stopping criterion based on Central Limit Theorem for multiple replications.

>>> k = Keister(Lattice(seed=7))
>>> sc = CubQMCCLT(k,abs_tol=.05)
>>> solution,data = sc.integrate()
>>> solution
array([1.38030146])
>>> data
MeanVarDataRep (AccumulateData Object)
    solution        1.380
    indv_error_bound 6.92e-04
    ci_low          1.380
    ci_high         1.381
    ci_comb_low     1.380
    ci_comb_high    1.381
    solution_comb   1.380
    flags_comb      0
    flags_indv      0
    n_total         2^(12)
    n               2^(8)
    replications    2^(4)
    time_integrate  ...
CubQMCCLT (StoppingCriterion Object)
    inflate         1.200
    alpha           0.010
    abs_tol         0.050
    rel_tol         0
    n_init          2^(8)
    n_max           2^(30)
Keister (Integrand Object)
Gaussian (TrueMeasure Object)
    mean            0
    covariance      2^(-1)
    decomp_type     PCA
Lattice (DiscreteDistribution Object)
    d               1
    dvec            0
    randomize       1
    order           natural
    entropy         7
    spawn_key       ()
>>> f = BoxIntegral(Lattice(3,seed=7), s=[-1,1])
>>> abs_tol = 1e-3
>>> sc = CubQMCCLT(f, abs_tol=abs_tol)
>>> solution,data = sc.integrate()
>>> solution
array([1.19023153, 0.96068581])
>>> data
MeanVarDataRep (AccumulateData Object)
    solution        [1.19  0.961]
    indv_error_bound [0.001 0.001]
    ci_low          [1.19 0.96]
    ci_high         [1.191 0.961]
    ci_comb_low     [1.19 0.96]
    ci_comb_high    [1.191 0.961]
    solution_comb   [1.19  0.961]
    flags_comb      [False False]
    flags_indv      [False False]
    n_total         2^(21)
    n               [131072.    512.]
    replications    2^(4)
    time_integrate  ...
CubQMCCLT (StoppingCriterion Object)
    inflate         1.200
    alpha           0.010
    abs_tol         0.001
    rel_tol         0
    n_init          2^(8)
    n_max           2^(30)
BoxIntegral (Integrand Object)
    s               [-1  1]
Uniform (TrueMeasure Object)
    lower_bound     0
    upper_bound     1
Lattice (DiscreteDistribution Object)
    d               3
    dvec            [0 1 2]
    randomize       1
    order           natural
    entropy         7
    spawn_key       ()
>>> sol3neg1 = -pi/4-1/2*log(2)+log(5+3*sqrt(3))
>>> sol31 = sqrt(3)/4+1/2*log(2+sqrt(3))-pi/24
>>> true_value = array([sol3neg1,sol31])
>>> (abs(true_value-solution)<abs_tol).all()
True
__init__(integrand, abs_tol=0.01, rel_tol=0.0, n_init=256.0, n_max=1073741824, inflate=1.2, alpha=0.01, replications=16.0, error_fun=<function CubQMCCLT.<lambda>>)
Parameters
  • integrand (Integrand) – an instance of Integrand

  • inflate (float) – inflation factor when estimating variance

  • alpha (float) – significance level for confidence interval

  • abs_tol (float) – absolute error tolerance

  • rel_tol (float) – relative error tolerance

  • n_max (int) – maximum number of samples

  • replications (int) – number of replications

  • error_fun – function taking in the approximate solution vector, absolute tolerance, and relative tolerance which returns the approximate error. Default indicates integration until either absolute OR relative tolerance is satisfied.

integrate()

See abstract method.

set_tolerance(abs_tol=None, rel_tol=None)

See abstract method.

Parameters
  • abs_tol (float) – absolute tolerance. Reset if supplied, ignored if not.

  • rel_tol (float) – relative tolerance. Reset if supplied, ignored if not.

Guaranteed MC Cubature

class qmcpy.stopping_criterion.cub_mc_g.CubMCG(integrand, abs_tol=0.01, rel_tol=0.0, n_init=1024.0, n_max=10000000000.0, inflate=1.2, alpha=0.01, control_variates=[], control_variate_means=[])

Stopping criterion with guaranteed accuracy.

>>> k = Keister(IIDStdUniform(2,seed=7))
>>> sc = CubMCG(k,abs_tol=.05)
>>> solution,data = sc.integrate()
>>> data
MeanVarData (AccumulateData Object)
    solution        1.807
    error_bound     0.050
    n_total         15256
    n               14232
    levels          1
    time_integrate  ...
CubMCG (StoppingCriterion Object)
    abs_tol         0.050
    rel_tol         0
    n_init          2^(10)
    n_max           10000000000
    inflate         1.200
    alpha           0.010
Keister (Integrand Object)
Gaussian (TrueMeasure Object)
    mean            0
    covariance      2^(-1)
    decomp_type     PCA
IIDStdUniform (DiscreteDistribution Object)
    d               2^(1)
    entropy         7
    spawn_key       ()
>>> dd = IIDStdUniform(1,seed=7)
>>> k = Keister(dd)
>>> cv1 = CustomFun(Uniform(dd),lambda x: sin(pi*x).sum(1))
>>> cv1mean = 2/pi
>>> cv2 = CustomFun(Uniform(dd),lambda x: (-3*(x-.5)**2+1).sum(1))
>>> cv2mean = 3/4
>>> sc1 = CubMCG(k,abs_tol=.05,control_variates=[cv1,cv2],control_variate_means=[cv1mean,cv2mean])
>>> sol,data = sc1.integrate()
>>> sol
1.384...

Original Implementation:

References

[1] Fred J. Hickernell, Lan Jiang, Yuewei Liu, and Art B. Owen, “Guaranteed conservative fixed width confidence intervals via Monte Carlo sampling,” Monte Carlo and Quasi-Monte Carlo Methods 2012 (J. Dick, F. Y. Kuo, G. W. Peters, and I. H. Sloan, eds.), pp. 105-128, Springer-Verlag, Berlin, 2014. DOI: 10.1007/978-3-642-41095-6_5

[2] Sou-Cheng T. Choi, Yuhan Ding, Fred J. Hickernell, Lan Jiang, Lluis Antoni Jimenez Rugama, Da Li, Jagadeeswaran Rathinavel, Xin Tong, Kan Zhang, Yizhi Zhang, and Xuan Zhou, GAIL: Guaranteed Automatic Integration Library (Version 2.3) [MATLAB Software], 2019. Available from http://gailgithub.github.io/GAIL_Dev/

Guarantee:

This algorithm attempts to calculate the mean, mu, of a random variable to a prescribed error tolerance, _tol_fun:= max(abstol,reltol*|mu|), with guaranteed confidence level 1-alpha. If the algorithm terminates without showing any warning messages and provides an answer tmu, then the follow inequality would be satisfied: \P(| mu - tmu | \\le tol\_fun) \\ge 1 - \\alpha.

__init__(integrand, abs_tol=0.01, rel_tol=0.0, n_init=1024.0, n_max=10000000000.0, inflate=1.2, alpha=0.01, control_variates=[], control_variate_means=[])
Parameters
  • integrand (Integrand) – an instance of Integrand

  • inflate – inflation factor when estimating variance

  • alpha – significance level for confidence interval

  • abs_tol – absolute error tolerance

  • rel_tol – relative error tolerance

  • n_init – initial number of samples

  • n_max – maximum number of samples

  • control_variates (list) – list of integrand objects to be used as control variates. Control variates are currently only compatible with single level problems. The same discrete distribution instance must be used for the integrand and each of the control variates.

  • control_variate_means (list) – list of means for each control variate

integrate()

See abstract method.

set_tolerance(abs_tol=None, rel_tol=None)

See abstract method.

Parameters
  • abs_tol (float) – absolute tolerance. Reset if supplied, ignored if not.

  • rel_tol (float) – relative tolerance. Reset if supplied, ignored if not.

CLT MC Cubature

class qmcpy.stopping_criterion.cub_mc_clt.CubMCCLT(integrand, abs_tol=0.01, rel_tol=0.0, n_init=1024.0, n_max=10000000000.0, inflate=1.2, alpha=0.01, control_variates=[], control_variate_means=[])

Stopping criterion based on the Central Limit Theorem.

>>> ao = AsianOption(IIDStdUniform(seed=7))
>>> sc = CubMCCLT(ao,abs_tol=.05)
>>> solution,data = sc.integrate()
>>> data
MeanVarData (AccumulateData Object)
    solution        1.519
    error_bound     0.046
    n_total         96028
    n               95004
    levels          1
    time_integrate  ...
CubMCCLT (StoppingCriterion Object)
    abs_tol         0.050
    rel_tol         0
    n_init          2^(10)
    n_max           10000000000
    inflate         1.200
    alpha           0.010
AsianOption (Integrand Object)
    volatility      2^(-1)
    call_put        call
    start_price     30
    strike_price    35
    interest_rate   0
    mean_type       arithmetic
    dim_frac        0
BrownianMotion (TrueMeasure Object)
    time_vec        1
    drift           0
    mean            0
    covariance      1
    decomp_type     PCA
IIDStdUniform (DiscreteDistribution Object)
    d               1
    entropy         7
    spawn_key       ()
>>> ao = AsianOption(IIDStdUniform(seed=7),multilevel_dims=[2,4,8])
>>> sc = CubMCCLT(ao,abs_tol=.05)
>>> solution,data = sc.integrate()
>>> dd = IIDStdUniform(1,seed=7)
>>> k = Keister(dd)
>>> cv1 = CustomFun(Uniform(dd),lambda x: sin(pi*x).sum(1))
>>> cv1mean = 2/pi
>>> cv2 = CustomFun(Uniform(dd),lambda x: (-3*(x-.5)**2+1).sum(1))
>>> cv2mean = 3/4
>>> sc1 = CubMCCLT(k,abs_tol=.05,control_variates=[cv1,cv2],control_variate_means=[cv1mean,cv2mean])
>>> sol,data = sc1.integrate()
>>> sol
1.381...
__init__(integrand, abs_tol=0.01, rel_tol=0.0, n_init=1024.0, n_max=10000000000.0, inflate=1.2, alpha=0.01, control_variates=[], control_variate_means=[])
Parameters
  • integrand (Integrand) – an instance of Integrand

  • inflate (float) – inflation factor when estimating variance

  • alpha (float) – significance level for confidence interval

  • abs_tol (float) – absolute error tolerance

  • rel_tol (float) – relative error tolerance

  • n_max (int) – maximum number of samples

  • control_variates (list) – list of integrand objects to be used as control variates. Control variates are currently only compatible with single level problems. The same discrete distribution instance must be used for the integrand and each of the control variates.

  • control_variate_means (list) – list of means for each control variate

integrate()

See abstract method.

set_tolerance(abs_tol=None, rel_tol=None)

See abstract method.

Parameters
  • abs_tol (float) – absolute tolerance. Reset if supplied, ignored if not.

  • rel_tol (float) – relative tolerance. Reset if supplied, ignored if not.

Continuation Multilevel QMC Cubature

class qmcpy.stopping_criterion.cub_qmc_ml_cont.CubQMCMLCont(integrand, abs_tol=0.05, alpha=0.01, rmse_tol=None, n_init=256.0, n_max=10000000000.0, replications=32.0, levels_min=2, levels_max=10, n_tols=10, tol_mult=1.6681005372000588, theta_init=0.5)

Stopping criterion based on continuation multi-level quasi-Monte Carlo.

>>> mlco = MLCallOptions(Lattice(seed=7))
>>> sc = CubQMCMLCont(mlco,abs_tol=.05)
>>> solution,data = sc.integrate()
>>> data
MLQMCData (AccumulateData Object)
    solution        10.421
    n_total         98304
    n_level         [2048.  256.  256.  256.  256.]
    levels          5
    mean_level      [10.054  0.183  0.102  0.054  0.028]
    var_level       [2.027e-04 5.129e-05 2.656e-05 1.064e-05 3.466e-06]
    bias_estimate   0.016
    time_integrate  ...
CubQMCMLCont (StoppingCriterion Object)
    rmse_tol        0.019
    n_init          2^(8)
    n_max           10000000000
    replications    2^(5)
    levels_min      2^(1)
    levels_max      10
    n_tols          10
    tol_mult        1.668
    theta_init      2^(-1)
    theta           2^(-3)
MLCallOptions (Integrand Object)
    option          european
    sigma           0.200
    k               100
    r               0.050
    t               1
    b               85
    level           0
Gaussian (TrueMeasure Object)
    mean            0
    covariance      1
    decomp_type     PCA
Lattice (DiscreteDistribution Object)
    d               1
    dvec            0
    randomize       1
    order           natural
    entropy         7
    spawn_key       ()

References

[1] https://github.com/PieterjanRobbe/MultilevelEstimators.jl

__init__(integrand, abs_tol=0.05, alpha=0.01, rmse_tol=None, n_init=256.0, n_max=10000000000.0, replications=32.0, levels_min=2, levels_max=10, n_tols=10, tol_mult=1.6681005372000588, theta_init=0.5)
Parameters
  • integrand (Integrand) – integrand with multi-level g method

  • abs_tol (float) – absolute tolerance

  • alpha (float) – uncertaintly level. If rmse_tol not supplied, then rmse_tol = abs_tol/norm.ppf(1-alpha/2)

  • rmse_tol (float) – root mean squared error If supplied (not None), then absolute tolerance and alpha are ignored in favor of the rmse tolerance

  • n_max (int) – maximum number of samples

  • replications (int) – number of replications on each level

  • levels_min (int) – minimum level of refinement >= 2

  • levels_max (int) – maximum level of refinement >= Lmin

  • n_tols (int) – number of coarser tolerances to run

  • tol_mult (float) – coarser tolerance multiplication factor

  • theta_init (float) – initial error splitting constant

integrate()

ABSTRACT METHOD to determine the number of samples needed to satisfy the tolerance.

Returns

tuple containing:
  • solution (float): approximation to the integral

  • data (AccumulateData): an AccumulateData object

Return type

tuple

set_tolerance(abs_tol=None, alpha=0.01, rmse_tol=None)

See abstract method.

Parameters
  • integrand (Integrand) – integrand with multi-level g method

  • abs_tol (float) – absolute tolerance. Reset if supplied, ignored if not.

  • alpha (float) – uncertaintly level. If rmse_tol not supplied, then rmse_tol = abs_tol/norm.ppf(1-alpha/2)

  • rel_tol (float) – relative tolerance. Reset if supplied, ignored if not. Takes priority over aboluste tolerance and alpha if supplied.

Multilevel QMC Cubature

class qmcpy.stopping_criterion.cub_qmc_ml.CubQMCML(integrand, abs_tol=0.05, alpha=0.01, rmse_tol=None, n_init=256.0, n_max=10000000000.0, replications=32.0, levels_min=2, levels_max=10)

Stopping criterion based on multi-level quasi-Monte Carlo.

>>> mlco = MLCallOptions(Lattice(seed=7))
>>> sc = CubQMCML(mlco,abs_tol=.05)
>>> solution,data = sc.integrate()
>>> data
MLQMCData (AccumulateData Object)
    solution        10.434
    n_total         172032
    n_level         [4096.  256.  256.  256.  256.  256.]
    levels          6
    mean_level      [10.053  0.183  0.102  0.054  0.028  0.014]
    var_level       [5.699e-05 5.129e-05 2.656e-05 1.064e-05 3.466e-06 1.113e-06]
    bias_estimate   0.007
    time_integrate  ...
CubQMCML (StoppingCriterion Object)
    rmse_tol        0.019
    n_init          2^(8)
    n_max           10000000000
    replications    2^(5)
MLCallOptions (Integrand Object)
    option          european
    sigma           0.200
    k               100
    r               0.050
    t               1
    b               85
    level           0
Gaussian (TrueMeasure Object)
    mean            0
    covariance      1
    decomp_type     PCA
Lattice (DiscreteDistribution Object)
    d               1
    dvec            0
    randomize       1
    order           natural
    entropy         7
    spawn_key       ()

References

[1] M.B. Giles and B.J. Waterhouse. ‘Multilevel quasi-Monte Carlo path simulation’. pp.165-181 in Advanced Financial Modelling, in Radon Series on Computational and Applied Mathematics, de Gruyter, 2009. http://people.maths.ox.ac.uk/~gilesm/files/radon.pdf

__init__(integrand, abs_tol=0.05, alpha=0.01, rmse_tol=None, n_init=256.0, n_max=10000000000.0, replications=32.0, levels_min=2, levels_max=10)
Parameters
  • integrand (Integrand) – integrand with multi-level g method

  • abs_tol (float) – absolute tolerance

  • alpha (float) – uncertaintly level. If rmse_tol not supplied, then rmse_tol = abs_tol/norm.ppf(1-alpha/2)

  • rmse_tol (float) – root mean squared error If supplied (not None), then absolute tolerance and alpha are ignored in favor of the rmse tolerance

  • n_max (int) – maximum number of samples

  • replications (int) – number of replications on each level

  • levels_min (int) – minimum level of refinement >= 2

  • levels_max (int) – maximum level of refinement >= Lmin

integrate()

See abstract method.

set_tolerance(abs_tol=None, alpha=0.01, rmse_tol=None)

See abstract method.

Parameters
  • integrand (Integrand) – integrand with multi-level g method

  • abs_tol (float) – absolute tolerance. Reset if supplied, ignored if not.

  • alpha (float) – uncertaintly level. If rmse_tol not supplied, then rmse_tol = abs_tol/norm.ppf(1-alpha/2)

  • rel_tol (float) – relative tolerance. Reset if supplied, ignored if not. Takes priority over aboluste tolerance and alpha if supplied.

Continuation Multilevel MC Cubature

class qmcpy.stopping_criterion.cub_mc_ml_cont.CubMCMLCont(integrand, abs_tol=0.05, alpha=0.01, rmse_tol=None, n_init=256.0, n_max=10000000000.0, levels_min=2, levels_max=10, n_tols=10, tol_mult=1.6681005372000588, theta_init=0.5)

Stopping criterion based on continuation multi-level monte carlo.

>>> mlco = MLCallOptions(IIDStdUniform(seed=7))
>>> sc = CubMCMLCont(mlco,abs_tol=.05)
>>> solution,data = sc.integrate()
>>> data
MLMCData (AccumulateData Object)
    solution        10.400
    n_total         1193331
    levels          2^(2)
    n_level         [1133772.   22940.    8676.    2850.]
    mean_level      [10.059  0.186  0.105  0.05 ]
    var_level       [1.959e+02 1.603e-01 4.567e-02 1.013e-02]
    cost_per_sample [1. 2. 4. 8.]
    alpha           0.942
    beta            1.992
    gamma           1.000
    time_integrate  ...
CubMCMLCont (StoppingCriterion Object)
    rmse_tol        0.019
    n_init          2^(8)
    levels_min      2^(1)
    levels_max      10
    n_tols          10
    tol_mult        1.668
    theta_init      2^(-1)
    theta           2^(-1)
MLCallOptions (Integrand Object)
    option          european
    sigma           0.200
    k               100
    r               0.050
    t               1
    b               85
    level           0
Gaussian (TrueMeasure Object)
    mean            0
    covariance      1
    decomp_type     PCA
IIDStdUniform (DiscreteDistribution Object)
    d               1
    entropy         7
    spawn_key       ()

References

[1] https://github.com/PieterjanRobbe/MultilevelEstimators.jl

__init__(integrand, abs_tol=0.05, alpha=0.01, rmse_tol=None, n_init=256.0, n_max=10000000000.0, levels_min=2, levels_max=10, n_tols=10, tol_mult=1.6681005372000588, theta_init=0.5)
Parameters
  • integrand (Integrand) – integrand with multi-level g method

  • abs_tol (float) – absolute tolerance. Reset if supplied, ignored if not.

  • alpha (float) – uncertaintly level. If rmse_tol not supplied, then rmse_tol = abs_tol/norm.ppf(1-alpha/2)

  • rel_tol (float) – relative tolerance. Reset if supplied, ignored if not. Takes priority over absolute tolerance and alpha if supplied.

  • n_init (int) – initial number of samples

  • n_max (int) – maximum number of samples

  • levels_min (int) – minimum level of refinement >= 2

  • levels_max (int) – maximum level of refinement >= Lmin

  • n_tols (int) – number of coarser tolerances to run

  • tol_mult (float) – coarser tolerance multiplication factor

  • theta_init (float) – initial error splitting constant

integrate()

ABSTRACT METHOD to determine the number of samples needed to satisfy the tolerance.

Returns

tuple containing:
  • solution (float): approximation to the integral

  • data (AccumulateData): an AccumulateData object

Return type

tuple

set_tolerance(abs_tol=None, alpha=0.01, rmse_tol=None)

See abstract method.

Parameters
  • abs_tol (float) – absolute tolerance. Reset if supplied, ignored if not.

  • alpha (float) – uncertaintly level. If rmse_tol not supplied, then rmse_tol = abs_tol/norm.ppf(1-alpha/2)

  • rel_tol (float) – relative tolerance. Reset if supplied, ignored if not. Takes priority over aboluste tolerance and alpha if supplied.

Multilevel MC Cubature

class qmcpy.stopping_criterion.cub_mc_ml.CubMCML(integrand, abs_tol=0.05, alpha=0.01, rmse_tol=None, n_init=256.0, n_max=10000000000.0, levels_min=2, levels_max=10, alpha0=- 1.0, beta0=- 1.0, gamma0=- 1.0)

Stopping criterion based on multi-level monte carlo.

>>> mlco = MLCallOptions(IIDStdUniform(seed=7))
>>> sc = CubMCML(mlco,abs_tol=.05)
>>> solution,data = sc.integrate()
>>> data
MLMCData (AccumulateData Object)
    solution        10.450
    n_total         1213658
    levels          7
    n_level         [1.173e+06 2.369e+04 1.174e+04 3.314e+03 1.144e+03 4.380e+02 1.690e+02]
    mean_level      [1.006e+01 1.856e-01 1.053e-01 5.127e-02 2.699e-02 1.558e-02 7.068e-03]
    var_level       [1.958e+02 1.596e-01 4.603e-02 1.057e-02 2.978e-03 8.701e-04 2.552e-04]
    cost_per_sample [ 1.  2.  4.  8. 16. 32. 64.]
    alpha           0.936
    beta            1.870
    gamma           1.000
    time_integrate  ...
CubMCML (StoppingCriterion Object)
    rmse_tol        0.019
    n_init          2^(8)
    levels_min      2^(1)
    levels_max      10
    theta           2^(-1)
MLCallOptions (Integrand Object)
    option          european
    sigma           0.200
    k               100
    r               0.050
    t               1
    b               85
    level           0
Gaussian (TrueMeasure Object)
    mean            0
    covariance      1
    decomp_type     PCA
IIDStdUniform (DiscreteDistribution Object)
    d               1
    entropy         7
    spawn_key       ()

Original Implementation:

References

[1] M.B. Giles. ‘Multi-level Monte Carlo path simulation’. Operations Research, 56(3):607-617, 2008. http://people.maths.ox.ac.uk/~gilesm/files/OPRE_2008.pdf.

__init__(integrand, abs_tol=0.05, alpha=0.01, rmse_tol=None, n_init=256.0, n_max=10000000000.0, levels_min=2, levels_max=10, alpha0=- 1.0, beta0=- 1.0, gamma0=- 1.0)
Parameters
  • integrand (Integrand) – integrand with multi-level g method

  • abs_tol (float) – absolute tolerance. Reset if supplied, ignored if not.

  • alpha (float) – uncertaintly level. If rmse_tol not supplied, then rmse_tol = abs_tol/norm.ppf(1-alpha/2)

  • rel_tol (float) – relative tolerance. Reset if supplied, ignored if not. Takes priority over absolute tolerance and alpha if supplied.

  • n_init (int) – initial number of samples

  • n_max (int) – maximum number of samples

  • levels_min (int) – minimum level of refinement >= 2

  • levels_max (int) – maximum level of refinement >= Lmin

  • alpha0 (float) – weak error is O(2^{-alpha0*level})

  • beta0 (float) – variance is O(2^{-bet0a*level})

  • gamma0 (float) – sample cost is O(2^{gamma0*level})

Note

if alpha, beta, gamma are not positive, then they will be estimated

integrate()

See abstract method.

set_tolerance(abs_tol=None, alpha=0.01, rmse_tol=None)

See abstract method.

Parameters
  • abs_tol (float) – absolute tolerance. Reset if supplied, ignored if not.

  • alpha (float) – uncertaintly level. If rmse_tol not supplied, then rmse_tol = abs_tol/norm.ppf(1-alpha/2)

  • rel_tol (float) – relative tolerance. Reset if supplied, ignored if not. Takes priority over aboluste tolerance and alpha if supplied.

Utilities

_images/util_err.png _images/util_warn.png
qmcpy.util.latnetbuilder_linker.latnetbuilder_linker(lnb_dir='./', out_dir='./', fout_prefix='lnb4qmcpy')
Parameters
  • lnb_dir (str) – relative path to directory where outputMachine.txt is stored e.g. ‘my_lnb/poly_lat/’

  • out_dir (str) – relative path to directory where output should be stored e.g. ‘my_lnb/poly_lat_qmcpy/’

  • fout_prefix (str) – start of output file name. e.g. ‘my_poly_lat_vec’

Returns

path to file which can be passed into QMCPy’s Lattice or Sobol’ in order to use

the linked latnetbuilder generating vector/matrix e.g. ‘my_poly_lat_vec.10.16.npy’

Return type

str

Adapted from latnetbuilder parser:

https://github.com/umontreal-simul/latnetbuilder/blob/master/python-wrapper/latnetbuilder/parse_output.py#L74