QMCPy Documentation

Discrete Distribution Class

Abstract Discrete Distribution Class

class qmcpy.discrete_distribution._discrete_distribution.DiscreteDistribution

Discrete Distribution abstract class. DO NOT INSTANTIATE.

gen_samples(*args)

ABSTRACT METHOD Generate samples from discrete distribution.

Parameters

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

Returns

n \times d array of samples, where n is number of observations and d is dimension

Return type

ndarray

plot(dim_x=0, dim_y=1, n=128, point_size=5, color='c', show=True, out=None)

Make a scatter plot from samples. Requires dimension >= 2.

Parameters

  • dim_x (int) – index of first dimension to be plotted on horizontal axis.

  • dim_y (int) – index of the second dimension to be plotted on vertical axis.

  • n (int) – number of samples to draw as self.gen_samples(n)

  • point_size (int) – ax.scatter(…,s=point_size)

  • color (str) – ax.scatter(…,color=color)

  • show (bool) – show plot or not?

  • out (str) – file name to output image. If None, the image is not output

Returns

fig,ax (tuple) from fig,ax = matplotlib.pyplot.subplots(…)

set_dimension(dimension)

ABSTRACT METHOD to reset the dimension of the problem.

Parameters

dimension (int) – new dimension to reset to

Note

May not be applicable to every discrete distribution (ex: CustomIIDDistribution).

set_seed(seed)

ABSTRACT METHOD to reset the seed of the problem.

Parameters

seed (int) – new seed for generator

Note

May not be applicable to every discrete distribution (ex: InverseCDFSampling)

Lattice

class qmcpy.discrete_distribution.lattice.lattice.Lattice(dimension=1, randomize=True, order='natural', seed=None, z_path=None)

Quasi-Random Lattice nets in base 2.

>>> l = Lattice(2,seed=7)
>>> l
Lattice (DiscreteDistribution Object)
    dimension       2^(1)
    randomize       1
    order           natural
    seed            7
    mimics          StdUniform
>>> l.gen_samples(4)
array([[0.076, 0.78 ],
       [0.576, 0.28 ],
       [0.326, 0.53 ],
       [0.826, 0.03 ]])
>>> l.set_dimension(3)
>>> l.gen_samples(n_min=4,n_max=8)
array([[0.563, 0.098, 0.353],
       [0.063, 0.598, 0.853],
       [0.813, 0.848, 0.103],
       [0.313, 0.348, 0.603]])
>>> 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, z_path=None)
Parameters

  • dimension (int) – dimension of samples

  • 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 (int) – seed the random number generator for reproducibility

  • z_path (str) – path to generating vector. z_path should be formatted like ‘lattice_vec.3600.20.npy’ where ‘name.d_max.m_max.npy’ and d_max is the maximum dimenion and 2^m_max is the max number samples supported

apply_randomization(x)

Apply a digital shift to the samples.

Parameters

x (ndarray) – un-randomized samples to be digitally shifted.

Returns

x with digital shift aplied.

Return type

ndarray

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

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_non_random (bool) – return both the samples with and without randomization

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.

set_dimension(dimension)

See abstract method.

set_seed(seed)

See abstract method.

Sobol’

class qmcpy.discrete_distribution.sobol.sobol.Sobol(dimension=1, randomize='LMS', graycode=False, seed=None, z_path=None, dim0=0)

Quasi-Random Sobol nets in base 2.

>>> s = Sobol(2,seed=7)
>>> s
Sobol (DiscreteDistribution Object)
    dimension       2^(1)
    randomize       1
    graycode        0
    seed            [61615 58564]
    mimics          StdUniform
    dim0            0
>>> s.gen_samples(4)
array([[0.783, 0.173],
       [0.128, 0.816],
       [0.72 , 0.664],
       [0.316, 0.334]])
>>> s.set_dimension(3)
>>> s.gen_samples(n_min=4,n_max=8)
array([[0.882, 0.932, 0.573],
       [0.035, 0.071, 0.379],
       [0.569, 0.418, 0.036],
       [0.474, 0.593, 0.982]])
>>> Sobol(dimension=2,randomize=False,graycode=True).gen_samples(n_min=2,n_max=4)
array([[0.75, 0.25],
       [0.25, 0.75]])
>>> Sobol(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/

[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', graycode=False, seed=None, z_path=None, dim0=0)
Parameters

  • dimension (int) – dimension of samples

  • randomize (bool) – If True, apply digital shift to generated samples. Note: Non-randomized Sobol’ sequence includes the origin.

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

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

  • z_path (str) – path to generating matricies. z_path sould be formatted like gen_mat.21201.32.msb.npy with name.d_max.m_max.msb_or_lsb.npy

  • dim0 (int) – first dimension

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

set_dim0(dim0)

Reset the first dimension

Parameters

dim0 (int) – first dimension

set_dimension(dimension)

Reset the dimension

Parameters

dimension (int) – new dimension

set_graycode(graycode)

Reset the graycode

Parameters

graycode (bool) – use graycode?

set_randomize(randomize)

Reset the randomization

Parameters

randomize (str) – randomization type. Either ‘LMS’: linear matrix scramble with digital shift ‘DS’: just the digital shift

set_seed(seeds)

Reset the seeds

Parameters

seeds (int/list/None) – new seeds

Halton

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

Quasi-Random Halton nets.

>>> h = Halton(2,seed=7)
>>> h.gen_samples(1)
array([[0.166, 0.363]])
>>> h.gen_samples(1)
array([[0.166, 0.363]])
>>> h.set_dimension(4)
>>> h.set_seed(8)
>>> h.gen_samples(2)
array([[0.323, 0.148, 0.623, 0.913],
       [0.823, 0.482, 0.223, 0.342]])
>>> h
Halton (DiscreteDistribution Object)
    dimension       2^(2)
    generalize      1
    randomize       1
    seed            2^(3)
    mimics          StdUniform

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, generalize=True, randomize=True, seed=None)
Parameters

  • dimension (int) – dimension of samples

  • generalize (bool) – generalize the Halton sequence?

  • randomize (bool/str) – If False, does not randomize Halton points. If True, will use ‘QRNG’ randomization as in [1]. Supports max dimension 360. You can also set radnomize=’QRNG’ or randomize=’Halton’ to explicitly select a randomization method. Halton OWEN supports up to 1000 dimensions.

  • seed (int) – seed the random number generator for reproducibility

Note

See References [1] and [2] for specific randomization methods and differences.

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

set_dimension(dimension)

See abstract method.

set_seed(seed)

See abstract method.

Korobov

class qmcpy.discrete_distribution.korobov.korobov.Korobov(dimension=1, generator=[1], randomize=True, seed=None, backend='QRNG')

Quasi-Random Korobov nets.

>>> k = Korobov(1,seed=7)
>>> k.gen_samples(2)
array([[0.982],
       [0.482]])
>>> k.gen_samples(2)
array([[0.982],
       [0.482]])
>>> k.set_dimension(3)
>>> k.set_seed(8)
>>> k.gen_samples(4)
array([[0.265, 0.153, 0.115],
       [0.515, 0.403, 0.365],
       [0.765, 0.653, 0.615],
       [0.015, 0.903, 0.865]])
>>> k
Korobov (DiscreteDistribution Object)
    dimension       3
    generator       1
    randomize       1
    seed            2^(3)
    mimics          StdUniform
    backend         QRNG
>>> Korobov(2,generator=[3,1]).gen_samples(4)
array([[0.807, 0.834],
       [0.557, 0.084],
       [0.307, 0.334],
       [0.057, 0.584]])

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.

__init__(dimension=1, generator=[1], randomize=True, seed=None, backend='QRNG')
Parameters

  • dimension (int) – dimension of samples

  • generator (ndarray of ints) – generator in {1,..,n-1} either a vector of length d or a single number (which is appropriately extended)

  • randomize (bool) – randomize the Korobov sequence? Note: Non-randomized Korobov sequence includes origin

  • seed (int) – seed the random number generator for reproducibility

  • backend (str) – backend generator must be “QRNG”

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

Generate samples

Parameters

n (int) – number of samples

Returns

n x d (dimension) array of samples

Return type

ndarray

set_dimension(dimension)

See abstract method.

set_seed(seed)

See abstract method.

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
IIDStdUniform (DiscreteDistribution Object)
    dimension       2^(1)
    seed            7
    mimics          StdUniform
>>> dd.gen_samples(4)
array([[0.076, 0.78 ],
       [0.438, 0.723],
       [0.978, 0.538],
       [0.501, 0.072]])
>>> dd.set_dimension(3)
>>> x = dd.gen_samples(5)
>>> x.shape
(5, 3)
__init__(dimension=1, seed=None)
Parameters

  • dimension (int) – dimension of samples

  • seed (int) – seed the random number generator for reproducibility

gen_samples(n)

Generate samples

Parameters

n (int) – Number of observations to generate

Returns

n x d (dimension) array of samples

Return type

ndarray

set_dimension(dimension)

See abstract class.

IID Standard Gaussian

class qmcpy.discrete_distribution.iid_std_gaussian.IIDStdGaussian(dimension=1, seed=None)

A wrapper around NumPy’s IID Standard Gaussian generator numpy.random.randn.

>>> dd = IIDStdGaussian(dimension=2,seed=7)
>>> dd
IIDStdGaussian (DiscreteDistribution Object)
    dimension       2^(1)
    seed            7
    mimics          StdGaussian
>>> dd.gen_samples(4)
array([[ 1.691e+00, -4.659e-01],
       [ 3.282e-02,  4.075e-01],
       [-7.889e-01,  2.066e-03],
       [-8.904e-04, -1.755e+00]])
>>> dd.set_dimension(3)
>>> x = dd.gen_samples(5)
>>> x.shape
(5, 3)
__init__(dimension=1, seed=None)
Parameters

  • dimension (int) – dimension of samples

  • seed (int) – seed the random number generator for reproducibility

gen_samples(n)

Generate samples

Parameters

n (int) – Number of observations to generate

Returns

n x d (dimension) array of samples

Return type

ndarray

set_dimension(dimension)

See abstract method.

Custom IID Distribution

class qmcpy.discrete_distribution.custom_iid_distribution.CustomIIDDistribution(custom_generator)

Custom IID Discrte Distribution.

>>> random.seed(7)
>>> cd = CustomIIDDistribution(lambda n: random.poisson(lam=5,size=(n,3)))
>>> cd
CustomIIDDistribution (DiscreteDistribution Object)
    dimension       None
>>> cd.gen_samples(2)
array([[6, 3, 3],
       [4, 6, 6]])
__init__(custom_generator)
Parameters

custom_generator (function) – custom generator of discrete distribution

gen_samples(n)

Generate samples

Parameters

n (int) – Number of observations to generate

Returns

n x d (dimension) array of samples

Return type

ndarray

Inverse CDF Sampling

class qmcpy.discrete_distribution.inverse_cdf_sampling.InverseCDFSampling(distribution_mimicking_uniform, inverse_cdf_fun=<function InverseCDFSampling.<lambda>>)

Sampling by inverse CDF transform applied to discrete distribution samples mimics standard uniform.

>>> _lambda = 1.5
>>> exp_pdf = lambda x,l=_lambda: l*exp(-l*x)
>>> exp_inverse_cdf = lambda u,l=_lambda: -log(1-u)/l
>>> exponential_measure = InverseCDFSampling(
...     distribution_mimicking_uniform = Sobol(dimension=2,seed=7),
...     inverse_cdf_fun = exp_inverse_cdf)
>>> exponential_measure
InverseCDFSampling (DiscreteDistribution Object)
    dimension       2^(1)
>>> exponential_measure.gen_samples(n_min=4,n_max=8)
array([[1.424, 1.792],
       [0.024, 0.049],
       [0.561, 0.361],
       [0.428, 0.599]])
Math for above example:

  • y \sim \text{Exp}(l)

  • \text{pdf }y: f(x) = l*\exp(-l*x)

  • \text{cdf }y: F(x)= 1-\exp(-l*x)

  • F^{-1}(u) = -\log(1-u)/l \sim \text{Exp}(l) \text{ for } u \sim \text{Uniform}(0,1)

__init__(distribution_mimicking_uniform, inverse_cdf_fun=<function InverseCDFSampling.<lambda>>)
Parameters

  • distribution_mimicking_uniform (DiscreteDistribution) – DiscreteDistribution instance which mimics standard uniform

  • inverse_cdf_fun (function) – function accepting samples mimicing uniform and applying inverse CDF transform. Must have 1 input argument accepting an ndarray of size n (observations) by d (diemsion)

gen_samples(*args, **kwargs)

Generate samples

Parameters

n (int) – Number of observations to generate

Returns

n x d (dimension) array of samples

Return type

ndarray

Acceptance Rejection Sampling

class qmcpy.discrete_distribution.acceptance_rejection_sampling.AcceptanceRejectionSampling(objective_pdf, measure_to_sample_from, c=None, c_est_draws=1024, c_est_inflation_factor=1.1, draws_multiple=inf)

Acceptance Rejection Sampling.

>>> def f(x):
...     x = x if x<.5 else 1-x
...     density = float(16*x)/3 if x<1./4 else 4./3
...     return density
>>> sampling_measure = Uniform(IIDStdUniform(1,seed=7))
>>> ars = AcceptanceRejectionSampling(
...     objective_pdf = f,
...     measure_to_sample_from = sampling_measure,
...     c_est_inflation_factor = 1)
>>> ars
AcceptanceRejectionSampling (DiscreteDistribution Object)
    dimension       1
    c               1.333
>>> x = ars.gen_samples(5)
>>> x
array([[0.575],
       [0.604],
       [0.144],
       [0.848],
       [0.129]])
>>> x.shape
(5, 1)
Define:

  • \rho(x) is pdf of measure we do not know how to generate from (mystery)

  • \nu(x) is the pdf of the continuous distribution which the discrete distribution mimics (known)

Prodecure:

  1. samples s_i from \nu(x)

  2. samples u_i from \text{Uniform}(0,1)

  3. if u_i \leq \rho(s_i)/(c*\nu(s_i)) then keep s_i

Note

If c is not supplied, then this algorithm conservitively estimates c by taking an initial samples of size n_0 combined with inflation factor w such that c \approx w*\max(\rho(s_i)/\nu(s_i) \text{ for } i=1,...,n_0)

__init__(objective_pdf, measure_to_sample_from, c=None, c_est_draws=1024, c_est_inflation_factor=1.1, draws_multiple=inf)
Parameters

  • objective_pdf (function) – pdf function of objective measure

  • measure_to_sample_from (TrueMeasure) – true measure we can sample from

  • c (float) – c such that for any x in the domain c*\nu(x) \geq \rho(x). If not supplied c will be estimated

  • c_est_draws (int) – initial number of draws used to estimate c

  • c_est_inflation_factor (float) – w in Note

  • draws_multiple (float) – will raise exception if drawing over n*draws_multiple samples when trying to get n samples

  • inflate_c_factor (float) – c = possibly inflate c to avoid underestimating

gen_samples(n)

Generate samples

Parameters

n (int) – Number of observations to generate

Returns

n x d (dimension) array of samples

Return type

ndarray

True Measure Class

Abstract Measure Class

class qmcpy.true_measure._true_measure.TrueMeasure

True Measure abstract class. DO NOT INSTANTIATE.

gen_samples(*args, **kwargs)

ABSTRACT METHOD to generate samples from the DiscreteDistribution object and transform them to mimic TrueMeasure samples.

Parameters

  • *args (tuple) – Ordered arguments to self.distribution.gen_samples

  • **kwrags (dict) – Keyword arguments to self.distribution.gen_samples

Returns

samples from the DiscreteDistribution transformed to mimic the TrueMeasure.

Return type

ndarray

Note

May not be applicable for all measures (ex: Lebesgue).

pdf(x)

Probability density function

Parameters

x (ndarray) – d (dimension) vector sample at which to evaluate the pdf

Note

May not be applicable for all measures (ex: Lebesgue).

plot(*args, **kwargs)

Create a plot relevant to the true measure object.

set_dimension(dimension)

ABSTRACT METHOD to reset the dimension of the problem. A wrapper around DiscreteDistribution.set_dimension.

Parameters

dimension (int) – new dimension to reset to

Note

May not be applicable for all measures (ex: Gaussian with covariance != k*eye(d) for scalar k)

Uniform

class qmcpy.true_measure.uniform.Uniform(distribution, lower_bound=0.0, upper_bound=1.0)
>>> u = Uniform(Sobol(2,seed=7),lower_bound=1,upper_bound=2)
>>> u
Uniform (TrueMeasure Object)
    lower_bound     [1 1]
    upper_bound     [2 2]
>>> u.gen_samples(n_min=4,n_max=8)
array([[1.882, 1.932],
       [1.035, 1.071],
       [1.569, 1.418],
       [1.474, 1.593]])
>>> u.set_dimension(4)
>>> u
Uniform (TrueMeasure Object)
    lower_bound     [1 1 1 1]
    upper_bound     [2 2 2 2]
>>> u.gen_samples(n_min=4,n_max=8)
array([[1.882, 1.932, 1.573, 1.07 ],
       [1.035, 1.071, 1.379, 1.6  ],
       [1.569, 1.418, 1.036, 1.889],
       [1.474, 1.593, 1.982, 1.422]])
>>> u2 = Uniform(Sobol(2),lower_bound=[-.5,0],upper_bound=[1,3])
>>> u2
Uniform (TrueMeasure Object)
    lower_bound     [-0.5  0. ]
    upper_bound     [1 3]
>>> u2.pdf([0,1])
0.2222222222222222
__init__(distribution, lower_bound=0.0, upper_bound=1.0)
Parameters

  • distribution (DiscreteDistribution) – DiscreteDistribution instance

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

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

gen_samples(*args, **kwargs)

See abstract method.

pdf(x)

See abstract class.

plot(dim_x=0, dim_y=1, n=128, point_size=5, color='c', show=True, out=None)

Make a scatter plot from samples. Requires dimension >= 2.

Parameters

  • dim_x (int) – index of first dimension to be plotted on horizontal axis.

  • dim_y (int) – index of the second dimension to be plotted on vertical axis.

  • n (int) – number of samples to draw as self.gen_samples(n)

  • point_size (int) – ax.scatter(…,s=point_size)

  • color (str) – ax.scatter(…,color=color)

  • show (bool) – show plot or not?

  • out (str) – file name to output image. If None, the image is not output

Returns

fig,ax from fig,ax = matplotlib.pyplot.subplots(…)

Return type

tuple

set_dimension(dimension)

See abstract method.

Gaussian

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

Normal Measure.

>>> dd = Sobol(2,seed=7)
>>> g = Gaussian(dd,mean=1,covariance=1./4)
>>> g
Gaussian (TrueMeasure Object)
    mean            1
    covariance      2^(-2)
    decomp_type     pca
>>> g.gen_samples(n_min=4,n_max=8)
array([[1.592, 1.745],
       [0.096, 0.267],
       [1.087, 0.897],
       [0.967, 1.117]])
>>> g.set_dimension(4)
>>> g.gen_samples(n_min=2,n_max=4)
array([[1.291, 1.211, 0.774, 1.483],
       [0.761, 0.785, 1.165, 0.826]])
>>> g2 = Gaussian(Sobol(2),mean=[1,2],covariance=[[1,.5],[.5,2]])
>>> g2
Gaussian (TrueMeasure Object)
    mean            [1 2]
    covariance      [[1.  0.5]
                    [0.5 2. ]]
    decomp_type     pca
>>> g2.pdf(array([0,0]))
array([[0.038]])
__init__(distribution, mean=0.0, covariance=1.0, decomp_type='PCA')
Parameters

  • distribution (DiscreteDistribution) – DiscreteDistribution instance

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

gen_samples(*args, **kwargs)

See abstract method.

pdf(x)

See abstract method.

plot(dim_x=0, dim_y=1, n=128, point_size=5, color='c', show=True, out=None)

Make a scatter plot from samples. Requires dimension >= 2.

Parameters

  • dim_x (int) – index of first dimension to be plotted on horizontal axis.

  • dim_y (int) – index of the second dimension to be plotted on vertical axis.

  • n (int) – number of samples to draw as self.gen_samples(n)

  • point_size (int) – ax.scatter(…,s=point_size)

  • color (str) – ax.scatter(…,color=color)

  • show (bool) – show plot or not?

  • out (str) – file name to output image. If None, the image is not output

Returns

fig,ax from fig,ax = matplotlib.pyplot.subplots(…)

Return type

tuple

set_dimension(dimension)

See abstract method.

Brownian Motion

class qmcpy.true_measure.brownian_motion.BrownianMotion(distribution, assembly_type='PCA', drift=0.0)

Geometric Brownian Motion.

>>> dd = Sobol(2,seed=7)
>>> bm = BrownianMotion(dd,drift=1)
>>> bm
BrownianMotion (TrueMeasure Object)
    time_vector     [0.5 1. ]
    drift           1
    assembly_type   pca
>>> bm.gen_samples(n_min=4,n_max=8)
array([[ 0.659,  2.495],
       [-0.043, -1.098],
       [ 0.681,  1.121],
       [ 0.374,  0.99 ]])
>>> bm.set_dimension(4)
>>> bm
BrownianMotion (TrueMeasure Object)
    time_vector     [0.25 0.5  0.75 1.  ]
    drift           1
    assembly_type   pca
__init__(distribution, assembly_type='PCA', drift=0.0)
Parameters

  • distribution (DiscreteDistribution) – DiscreteDistribution instance

  • assembly_type (str) – assembly type type either “Diff” for time differencing, “PCA” for principal component analysis, or “Bridge” for Brownian Bridge.

  • drift (float) – mean shift for importance sampling.

gen_samples(*args, **kwargs)

See abstract method.

plot(n=32, show=True, out=None)

Plot Brownian Motion value vs time

Parameters

  • n (int) – self.gen_samples(n)

  • show (bool) – show the plot?

  • out (str) – file name to output image. If None, the image is not output

Returns

fig,ax from fig,ax = pyplot.subplots()

Return type

tuple

set_dimension(dimension)

See abstract method.

Note

Monitoring times are evenly spaced as linspace(1/dimension,1,dimension)

Lebesgue

class qmcpy.true_measure.lebesgue.Lebesgue(distribution, lower_bound=0.0, upper_bound=1.0)
>>> dd = Sobol(2,seed=7)
>>> Lebesgue(dd,lower_bound=[-1,0],upper_bound=[1,3])
Lebesgue (TrueMeasure Object)
    lower_bound     [-1  0]
    upper_bound     [1 3]
>>> Lebesgue(dd,lower_bound=-inf,upper_bound=inf)
Lebesgue (TrueMeasure Object)
    lower_bound     [-inf -inf]
    upper_bound     [inf inf]
__init__(distribution, lower_bound=0.0, upper_bound=1.0)
Parameters

  • distribution (DiscreteDistribution) – DiscreteDistribution instance

  • lower_bound (float or inf) – lower bound of integration

  • upper_bound (float or inf) – upper bound of integration

Importance Sampling

class qmcpy.true_measure.importance_sampling.ImportanceSampling(objective_pdf, measure_to_sample_from)

Importance Sampling.

>>> def quarter_circle_uniform_pdf(x):
...     x1,x2 = x
...     if sqrt(x1**2+x2**2)<1 and x1>=0 and x2>=0:
...         return 4./pi # 1./(pi*(1**2)/4)
...     else:
...         return 0. # outside of quarter circle
>>> tm = ImportanceSampling(
...     objective_pdf = quarter_circle_uniform_pdf,
...     measure_to_sample_from = Uniform(Lattice(2,seed=7)))
__init__(objective_pdf, measure_to_sample_from)
Parameters

  • objective_pdf (function) – pdf function of objective measure

  • measure_to_sample_from (TrueMeasure) – true measure we can sample from

Identical to what Discrete Distribution Mimics

class qmcpy.true_measure.identical_to_discrete.IdentitalToDiscrete(distribution)

For when Discrete Distribution samples already mimic the TrueMeasure. AKA: when g, the original integrand, and f, the transformed integrand are such that: g(x) = f(x) for x ~ DiscreteDistribution.

>>> dd = Sobol(2,seed=7,randomize=False,graycode=False)
>>> itd = IdentitalToDiscrete(dd)
>>> itd.gen_samples(4)
array([[0.  , 0.  ],
       [0.5 , 0.5 ],
       [0.25, 0.75],
       [0.75, 0.25]])
__init__(distribution)
Parameters

distribution (DiscreteDistribution) – DiscreteDistribution instance

gen_samples(*args, **kwargs)

See abstract method.

plot(*args, **kwargs)

Throw error if trying to plot.

Integrand Class

Abstract Integrand Class

class qmcpy.integrand._integrand.Integrand

Integrand abstract class. DO NOT INSTANTIATE.

g(x)

ABSTRACT METHOD for original integrand to be integrated.

Parameters

  • x (ndarray) – n samples by d dimension array of samples generated according to the true measure.

  • l (int) – OPTIONAL input for multi-level integrands. The level to generate at. Note that the dimension of x is determined by the _dim_at_level method for multi-level methods.

Returns

n vector of function evaluations

Return type

ndarray

period_transform(ptransform)

Computes the periodization transform for the given function values

plot(*args, **kwargs)

Create a plot relevant to the true measure object.

Keister Function

class qmcpy.integrand.keister.Keister(measure)

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

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

>>> dd = Sobol(2,seed=7)
>>> m = Gaussian(dd,covariance=1./2)
>>> k = Keister(m)
>>> x = dd.gen_samples(2**10)
>>> y = k.f(x)
>>> y.mean()
1.80...

References

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

__init__(measure)
Parameters

measure (TrueMeasure) – a TrueMeasure instance

g(x)

See abstract method.

plot(projection_dims=[0], n=128, point_size=5, color='c', show=True, out=None)

Make a scatter plot from samples. Requires dimension >= 2.

Parameters

  • projection_dims (list of ints) – dimensions to project onto individual dimensions. For example: projection_dims=[0,1] will make 2 plots. One with y vs x_0 and one with y vs x_1.

  • n (int) – number of samples to draw as self.gen_samples(n)

  • point_size (int) – ax.scatter(…,s=point_size)

  • color (str) – ax.scatter(…,color=color)

  • show (bool) – show plot or not?

  • out (str) – file name to output image. If None, the image is not output

Returns

fig,ax from fig,ax = matplotlib.pyplot.subplots(…)

Return type

tuple

Custom Function

class qmcpy.integrand.custom_fun.CustomFun(measure, custom_fun)

Custom user-supplied function handle.

>>> dd = Sobol(2,seed=7)
>>> m = Gaussian(dd,mean=[1,2])
>>> cf = CustomFun(m,lambda x: x[:,0]**2*x[:,1])
>>> x = dd.gen_samples(2**10)
>>> y = cf.f(x)
>>> y.mean()
4.00...
__init__(measure, custom_fun)
Parameters

  • measure (TrueMeasure) – a TrueMeasure instance

  • custom_fun (function) – a function evaluating samples (nxd) -> (nx1). See g method.

g(x)

See abstract method.

European Option

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

European financial option.

>>> dd = Sobol(4,seed=17)
>>> m = BrownianMotion(dd,drift=-1)
>>> eo = EuropeanOption(m,call_put='put')
>>> eo
EuropeanOption (Integrand Object)
    volatility      2^(-1)
    call_put        put
    start_price     30
    strike_price    35
    interest_rate   0
>>> x = dd.gen_samples(2**12)
>>> y = eo.f(x)
>>> y.mean()
9.2...
__init__(measure, volatility=0.5, start_price=30, strike_price=35, interest_rate=0, call_put='call')
Parameters

  • measure (TrueMeasure) – A BrownianMotion TrueMeasure object

  • 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

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

g(x)

See abstract method.

get_exact_value()

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

Returns

fair price

Return type

float

plot(n=32, show=True, out=None)

Plot European Option price vs time.

Parameters

  • n (int) – self.gen_samples(n)

  • show (bool) – show the plot?

  • out (str) – file name to output image. If None, the image is not output

Returns

fig,ax from fig,ax = pyplot.subplots()

Return type

tuple

Asian Option

class qmcpy.integrand.asian_option.AsianOption(measure, volatility=0.5, start_price=30.0, strike_price=35.0, interest_rate=0.0, call_put='call', mean_type='arithmetic', multi_level_dimensions=None)

Asian financial option.

>>> dd = Sobol(4,seed=7)
>>> m = BrownianMotion(dd)
>>> ac = AsianOption(m)
>>> ac
AsianOption (Integrand Object)
    volatility      2^(-1)
    call_put        call
    start_price     30
    strike_price    35
    interest_rate   0
    mean_type       arithmetic
    dimensions      2^(2)
    dim_fracs       0
>>> x = dd.gen_samples(2**10)
>>> y = ac.f(x)
>>> y.mean()
1.77...
>>> dd2 = Sobol(seed=7)
>>> m2 = BrownianMotion(dd2,drift=1)
>>> level_dims = [2,4,8]
>>> ac2 = AsianOption(m2,multi_level_dimensions=level_dims)
>>> ac2
AsianOption (Integrand Object)
    volatility      2^(-1)
    call_put        call
    start_price     30
    strike_price    35
    interest_rate   0
    mean_type       arithmetic
    dimensions      [2 4 8]
    dim_fracs       [0. 2. 2.]
>>> y2 = 0
>>> for l in range(len(level_dims)):
...     new_dim = ac2._dim_at_level(l)
...     m2.set_dimension(new_dim)
...     x2 = dd2.gen_samples(2**10)
...     y2 += ac2.f(x2,l=l).mean()
>>> y2
1.78...
__init__(measure, volatility=0.5, start_price=30.0, strike_price=35.0, interest_rate=0.0, call_put='call', mean_type='arithmetic', multi_level_dimensions=None)
Parameters

  • measure (TrueMeasure) – A BrownianMotion TrueMeasure object

  • 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

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

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

g(x, l=0)

See abstract method.

plot(n=32, show=True, out=None)

Plot Asian Option price vs time. Does not work for multi-level Asian option.

Parameters

  • n (int) – self.gen_samples(n)

  • show (bool) – show the plot?

  • out (str) – file name to output image. If None, the image is not output

Returns

fig,ax from fig,ax = pyplot.subplots()

Return type

tuple

Multilevel Call Options with Milstein Discretization

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

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

>>> dd = Sobol(seed=7)
>>> m = Gaussian(dd)
>>> mlco = MLCallOptions(m)
>>> mlco
MLCallOptions (Integrand Object)
    option          european
    sigma           0.200
    k               100
    r               0.050
    t               1
    b               85
>>> y = 0
>>> for level in range(4):
...     new_dim = mlco._dim_at_level(level)
...     m.set_dimension(new_dim)
...     x = dd.gen_samples(2**10)
...     sums,cost = mlco.f(x,l=level)
...     y += sums[0]/2**10
>>> y
10.39...

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__(measure, option='european', volatility=0.2, start_strike_price=100.0, interest_rate=0.05, t_final=1.0)
Parameters

  • measure (TrueMeasure) – A BrownianMotion TrueMeasure object

  • 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

g(samples, l)
Parameters

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

  • l (int) – level

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.linear.Linear(measure)

f(\boldsymbol{x}) = \sum_{i=1}^d x_i

>>> dd = Sobol(100,seed=7)
>>> m = Gaussian(dd,mean=(-1)**arange(100),covariance=1./3)
>>> l = Linear(m)
>>> x = dd.gen_samples(2**10)
>>> y = l.f(x)
>>> y.mean()
-0.002...
__init__(measure)
Parameters

measure (TrueMeasure) – a TrueMeasure instance

g(x)

See abstract method.

plot(projection_dims=[0], n=128, point_size=5, color='c', show=True, out=None)

Make a scatter plot from samples. Requires dimension >= 2.

Parameters

  • projection_dims (list of ints) – dimensions to project onto individual dimensions. For example: projection_dims=[0,1] will make 2 plots. One with y vs x_0 and one with y vs x_1.

  • n (int) – number of samples to draw as self.gen_samples(n)

  • point_size (int) – ax.scatter(…,s=point_size)

  • color (str) – ax.scatter(…,color=color)

  • show (bool) – show plot or not?

  • out (str) – file name to output image. If None, the image is not output

Returns

fig,ax from fig,ax = matplotlib.pyplot.subplots(…)

Return type

tuple

Stopping Criterion Algorithms

Abstract Stopping Criterion Class

class qmcpy.stopping_criterion._stopping_criterion.StoppingCriterion(distribution, integrand, allowed_levels, allowed_distribs)

Stopping Criterion abstract class. DO NOT INSTANTIATE.

__init__(distribution, integrand, allowed_levels, allowed_distribs)
Parameters

  • distribution (DiscreteDistribution) – a DiscreteDistribution

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

  • allowed_distribs (list) – list of names (strings) of compatible distributions

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

plot(*args, **kwargs)

Create a plot relevant to the stopping criterion object.

set_tolerance(*args, **kwargs)

ABSTRACT METHOD to reset the absolute tolerance.

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)

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(Gaussian(Lattice(2,seed=7),covariance=1./2))
>>> sc = CubQMCLatticeG(k,abs_tol=.05)
>>> solution,data = sc.integrate()
>>> solution
1.8079801428928022
>>> data
Solution: 1.8080
Keister (Integrand Object)
Lattice (DiscreteDistribution Object)
    dimension       2^(1)
    randomize       1
    order           natural
    seed            7
    mimics          StdUniform
Gaussian (TrueMeasure Object)
    mean            0
    covariance      2^(-1)
    decomp_type     pca
CubQMCLatticeG (StoppingCriterion Object)
    abs_tol         0.050
    rel_tol         0
    n_init          2^(10)
    n_max           2^(35)
LDTransformData (AccumulateData Object)
    n_total         2^(10)
    solution        1.808
    error_bound     0.005
    time_integrate  ...

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

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 Sobol Cubature (qMC)

class qmcpy.stopping_criterion.cub_qmc_sobol_g.CubQMCSobolG(integrand, abs_tol=0.01, rel_tol=0.0, n_init=1024.0, n_max=34359738368.0, fudge=<function CubQMCSobolG.<lambda>>, check_cone=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(Gaussian(Sobol(2,seed=7),covariance=1./2))
>>> sc = CubQMCSobolG(k,abs_tol=.05)
>>> solution,data = sc.integrate()
>>> solution
1.80...
>>> data
Solution: 1.8068
Keister (Integrand Object)
Sobol (DiscreteDistribution Object)
    dimension       2^(1)
    randomize       1
    graycode        0
    seed            [61615 58564]
    mimics          StdUniform
    dim0            0
Gaussian (TrueMeasure Object)
    mean            0
    covariance      2^(-1)
    decomp_type     pca
CubQMCSobolG (StoppingCriterion Object)
    abs_tol         0.050
    rel_tol         0
    n_init          2^(10)
    n_max           2^(35)
LDTransformData (AccumulateData Object)
    n_total         2^(10)
    solution        1.807
    error_bound     0.007
    time_integrate  ...

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 CubQMCSobolG.<lambda>>, check_cone=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

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.

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)

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

>>> mlco = MLCallOptions(Gaussian(Lattice(seed=7)))
>>> sc = CubQMCML(mlco,abs_tol=.05)
>>> solution,data = sc.integrate()
>>> solution
10.446333077585034
>>> data
Solution: 10.4463
MLCallOptions (Integrand Object)
    option          european
    sigma           0.200
    k               100
    r               0.050
    t               1
    b               85
Lattice (DiscreteDistribution Object)
    dimension       2^(6)
    randomize       1
    order           natural
    seed            854306
    mimics          StdUniform
Gaussian (TrueMeasure Object)
    mean            0
    covariance      1
    decomp_type     pca
CubQMCML (StoppingCriterion Object)
    rmse_tol        0.019
    n_init          2^(8)
    n_max           10000000000
    replications    2^(5)
MLQMCData (AccumulateData Object)
    levels          7
    dimensions      [ 1.  2.  4.  8. 16. 32. 64.]
    n_level         [4096.  512.  256.  256.  256.  256.  256.]
    mean_level      [1.006e+01 1.852e-01 1.040e-01 5.332e-02 2.751e-02 1.399e-02 6.994e-03]
    var_level       [6.143e-05 3.535e-05 2.896e-05 1.767e-05 2.898e-06 9.919e-07 3.260e-07]
    bias_estimate   0.007
    n_total         188416
    time_integrate  ...

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

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.

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)

Stopping criterion based on Central Limit Theorem for multiple replications.

>>> k = Keister(Gaussian(Lattice(seed=7),covariance=1./2))
>>> sc = CubQMCCLT(k,abs_tol=.05)
>>> solution,data = sc.integrate()
>>> solution
1.3798619783658828
>>> data
Solution: 1.3799
Keister (Integrand Object)
Lattice (DiscreteDistribution Object)
    dimension       1
    randomize       1
    order           natural
    seed            1092
    mimics          StdUniform
Gaussian (TrueMeasure Object)
    mean            0
    covariance      2^(-1)
    decomp_type     pca
CubQMCCLT (StoppingCriterion Object)
    inflate         1.200
    alpha           0.010
    abs_tol         0.050
    rel_tol         0
    n_init          2^(8)
    n_max           2^(30)
MeanVarDataRep (AccumulateData Object)
    replications    2^(4)
    solution        1.380
    sighat          0.001
    n_total         2^(12)
    error_bound     8.30e-04
    confid_int      [1.379 1.381]
    time_integrate  ...
__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)
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

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.

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.0, levels_max=10.0, alpha0=- 1.0, beta0=- 1.0, gamma0=- 1.0)

Stopping criterion based on multi-level monte carlo.

>>> mlco = MLCallOptions(Gaussian(IIDStdGaussian(seed=7)))
>>> sc = CubMCML(mlco,abs_tol=.05)
>>> solution,data = sc.integrate()
>>> solution
10.445441712933263
>>> data
Solution: 10.4454
MLCallOptions (Integrand Object)
    option          european
    sigma           0.200
    k               100
    r               0.050
    t               1
    b               85
IIDStdGaussian (DiscreteDistribution Object)
    dimension       2^(6)
    seed            7
    mimics          StdGaussian
Gaussian (TrueMeasure Object)
    mean            0
    covariance      1
    decomp_type     pca
CubMCML (StoppingCriterion Object)
    rmse_tol        0.019
    n_init          2^(8)
    levels_min      2^(1)
    levels_max      10
    theta           2^(-2)
MLMCData (AccumulateData Object)
    levels          7
    dimensions      [ 1.  2.  4.  8. 16. 32. 64.]
    n_level         [7.804e+05 1.533e+04 6.633e+03 2.077e+03 7.560e+02 2.730e+02 9.600e+01]
    mean_level      [1.006e+01 1.848e-01 1.014e-01 5.138e-02 2.472e-02 1.452e-02 7.657e-03]
    var_level       [1.963e+02 1.515e-01 4.124e-02 1.109e-02 2.901e-03 6.799e-04 1.848e-04]
    cost_per_sample [ 1.  2.  4.  8. 16. 32. 64.]
    n_total         805984
    alpha           0.927
    beta            1.946
    gamma           1.000
    time_integrate  ...

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.0, levels_max=10.0, 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 aboluste 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.

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)

Stopping criterion with guaranteed accuracy.

>>> k = Keister(Gaussian(IIDStdUniform(2,seed=7),covariance=1./2))
>>> sc = CubMCG(k,abs_tol=.05)
>>> solution,data = sc.integrate()
>>> solution
1.803926962264685
>>> data
Solution: 1.8039
Keister (Integrand Object)
IIDStdUniform (DiscreteDistribution Object)
    dimension       2^(1)
    seed            7
    mimics          StdUniform
Gaussian (TrueMeasure Object)
    mean            0
    covariance      2^(-1)
    decomp_type     pca
CubMCG (StoppingCriterion Object)
    inflate         1.200
    alpha           0.010
    abs_tol         0.050
    rel_tol         0
    n_init          2^(10)
    n_max           10000000000
MeanVarData (AccumulateData Object)
    levels          1
    solution        1.804
    n               13473
    n_total         14497
    error_bound     0.050
    confid_int      [1.754 1.854]
    time_integrate  ...

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

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)

Stopping criterion based on the Central Limit Theorem.

>>> k = Keister(Gaussian(IIDStdGaussian(2,seed=7),covariance=1./2))
>>> sc = CubMCCLT(k,abs_tol=.05)
>>> solution,data = sc.integrate()
>>> solution
1.834674149189029
>>> data
Solution: 1.8347
Keister (Integrand Object)
IIDStdGaussian (DiscreteDistribution Object)
    dimension       2^(1)
    seed            7
    mimics          StdGaussian
Gaussian (TrueMeasure Object)
    mean            0
    covariance      2^(-1)
    decomp_type     pca
CubMCCLT (StoppingCriterion Object)
    inflate         1.200
    alpha           0.010
    abs_tol         0.050
    rel_tol         0
    n_init          2^(10)
    n_max           10000000000
MeanVarData (AccumulateData Object)
    levels          1
    solution        1.835
    n               5826
    n_total         6850
    error_bound     0.050
    confid_int      [1.785 1.885]
    time_integrate  ...
>>> ac = AsianOption(
...     measure = BrownianMotion(IIDStdGaussian()),
...     multi_level_dimensions = [2,4,8])
>>> sc = CubMCCLT(ac,abs_tol=.05)
>>> solution,data = sc.integrate()
__init__(integrand, abs_tol=0.01, rel_tol=0.0, n_init=1024.0, n_max=10000000000.0, inflate=1.2, alpha=0.01)
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

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.

Accumulate Data Class

Abstract Accumulate Data Class

class qmcpy.accumulate_data._accumulate_data.AccumulateData

Accumulated Data abstract class. DO NOT INSTANTIATE.

__init__()

Initialize data instance

update_data()

ABSTRACT METHOD to update the accumulated data.

LD Sequence Transform Data (qMC)

class qmcpy.accumulate_data.ld_transform_data.LDTransformData(stopping_criterion, integrand, basis_transform, m_min, m_max, fudge, check_cone)

Update and store transformation data based on low-discrepancy sequences. See the stopping criterion that utilize this object for references.

__init__(stopping_criterion, integrand, basis_transform, m_min, m_max, fudge, check_cone)
Parameters

  • stopping_criterion (StoppingCriterion) – a StoppingCriterion instance

  • integrand (Integrand) – an Integrand instance

  • basis_transform (method) – Transform ynext, combine with y, and then transform all points. For cub_lattice this is Fast Fourier Transform (FFT). For cub_sobol this is Fast Walsh Transform (FWT)

  • m_min (int) – initial n == 2^m_min

  • m_max (int) – max n == 2^m_max

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

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

update_data()

See abstract method.

Mean Variance qMC Data (for Replications)

class qmcpy.accumulate_data.mean_var_data_rep.MeanVarDataRep(stopping_criterion, integrand, n_init, replications)

Update and store mean and variance estimates with repliations. See the stopping criterion that utilize this object for references.

__init__(stopping_criterion, integrand, n_init, replications)
Parameters

  • stopping_criterion (StoppingCriterion) – a StoppingCriterion instance

  • integrand (Integrand) – an Integrand instance

  • n_init (int) – initial number of samples

  • replications (int) – number of replications

update_data()

See abstract method.

Multilevel qMC Data

class qmcpy.accumulate_data.mlqmc_data.MLQMCData(stopping_criterion, integrand, n_init, replications)

Accumulated data for quasi-random sequence calculations, and store multi-level, multi-replications mean, variance, and cost values. See the stopping criterion that utilize this object for references.

__init__(stopping_criterion, integrand, n_init, replications)

Initialize data instance

Parameters

  • stopping_criterion (StoppingCriterion) – a StoppingCriterion instance

  • integrand (Integrand) – an Integrand instance

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

update_data()

See abstract method.

Multilevel MC Data

class qmcpy.accumulate_data.mlmc_data.MLMCData(stopping_criterion, integrand, levels_init, n_init, alpha0, beta0, gamma0)

Accumulated data for IIDDistribution calculations, and store multi-level mean, variance, and cost values. See the stopping criterion that utilize this object for references.

__init__(stopping_criterion, integrand, levels_init, n_init, alpha0, beta0, gamma0)

Initialize data instance

Parameters

  • stopping_criterion (StoppingCriterion) – a StoppingCriterion instance

  • integrand (Integrand) – an Integrand instance

  • levels_init (int) – initial number of levels

  • n_init (int) – initial number of samples per level

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

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

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

update_data()

See abstract method.

Mean Variance MC Data

class qmcpy.accumulate_data.mean_var_data.MeanVarData(stopping_criterion, integrand, n_init)

Update and store mean and variance estimates.

__init__(stopping_criterion, integrand, n_init)
Parameters

  • stopping_criterion (StoppingCriterion) – a StoppingCriterion instance

  • integrand (Integrand) – an Integrand instance

  • n_init (int) – initial number of samples

update_data()

See abstract method.