Skip to content

Stopping Criteria

UML Overview

AbstractStoppingCriterion

Bases: object

Parameters:

Name Type Description Default
allowed_distribs list

Allowed discrete distribution classes.

 required
allow_vectorized_integrals bool

Whether or not to allow integrands with vectorized outputs, i.e., those with integrand.d_indv!=().

 required
Source code in qmcpy/stopping_criterion/abstract_stopping_criterion.py 
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
def __init__(self, allowed_distribs, allow_vectorized_integrals):
    """
    Args:
        allowed_distribs (list): Allowed discrete distribution classes.
        allow_vectorized_integrals (bool): Whether or not to allow integrands with vectorized outputs, 
            i.e., those with `integrand.d_indv!=()`. 
    """
    sname = type(self).__name__
    prefix = 'A concrete implementation of StoppingCriterion must have '
    # integrand check
    if (not hasattr(self, 'integrand')) or (not isinstance(self.integrand,AbstractIntegrand)):
        raise ParameterError(prefix + 'self.integrand, an Integrand instance')
    # true measure check
    if (not hasattr(self, 'true_measure')) or (self.true_measure!=self.integrand.true_measure):
        raise ParameterError(prefix + 'self.true_measure=self.integrand.true_measure')
    # discrete distribution check
    if (not hasattr(self, 'discrete_distrib')) or (self.discrete_distrib!=self.integrand.discrete_distrib):
        raise ParameterError(prefix + 'self.discrete_distrib=self.integrand.discrete_distrib')
    if not isinstance(self.discrete_distrib,tuple(allowed_distribs)):
        raise DistributionCompatibilityError('%s must have a DiscreteDistribution in %s'%(sname,str(allowed_distribs)))
    if (not allow_vectorized_integrals) and len(self.integrand.d_indv)>0:
        raise ParameterError('Vectorized integrals (with d_indv!=() outputs per sample) are not supported by this stopping criterion')
    # parameter checks
    if not hasattr(self,'parameters'):
        self.parameters = []

integrate 

integrate()

Abstract method to determine the number of samples needed to satisfy the tolerance(s).

Returns:

Name Type Description
solution Union[float, ndarray]

Approximation to the integral with shape integrand.d_comb.
data Data

A data object.
Source code in qmcpy/stopping_criterion/abstract_stopping_criterion.py 
34
35
36
37
38
39
40
41
42
def integrate(self):
    """
    *Abstract method* to determine the number of samples needed to satisfy the tolerance(s).

    Returns:
        solution (Union[float,np.ndarray]): Approximation to the integral with shape `integrand.d_comb`.
        data (Data): A data object.
    """
    raise MethodImplementationError(self, 'integrate')

set_tolerance 

set_tolerance(abs_tol=None, rel_tol=None, rmse_tol=None)

Reset the tolerances.

Parameters:

Name Type Description Default
abs_tol float

Absolute tolerance (if applicable). Reset if supplied, ignored otherwise.

None
rel_tol float

Relative tolerance (if applicable). Reset if supplied, ignored otherwise.

None
rmse_tol float

RMSE tolerance (if applicable). Reset if supplied, ignored if not. If rmse_tol is not supplied but abs_tol is, then rmse_tol = abs_tol / norm.ppf(1-alpha/2).

 None
Source code in qmcpy/stopping_criterion/abstract_stopping_criterion.py 
59
60
61
62
63
64
65
66
67
68
69
def set_tolerance(self, abs_tol=None, rel_tol=None, rmse_tol=None):
    """
    Reset the tolerances.

    Args:
        abs_tol (float): Absolute tolerance (if applicable). Reset if supplied, ignored otherwise. 
        rel_tol (float): Relative tolerance (if applicable). Reset if supplied, ignored otherwise. 
        rmse_tol (float): RMSE tolerance (if applicable). Reset if supplied, ignored if not. 
            If `rmse_tol` is not supplied but `abs_tol` is, then `rmse_tol = abs_tol / norm.ppf(1-alpha/2)`. 
    """
    raise MethodImplementationError(self, 'integrate')

QMC

CubQMCNetG

Bases: AbstractCubQMCLDG

Quasi-Monte Carlo stopping criterion using digital net cubature with guarantees for cones of functions with a predictable decay in the Walsh coefficients.

Examples:

>>> k = Keister(DigitalNetB2(seed=7))
>>> sc = CubQMCNetG(k,abs_tol=1e-3,rel_tol=0,check_cone=True)
>>> solution,data = sc.integrate()
>>> solution
array(1.38046669)
>>> data
Data (Data)
    solution        1.380
    comb_bound_low  1.380
    comb_bound_high 1.381
    comb_bound_diff 6.72e-04
    comb_flags      1
    n_total         2^(11)
    n               2^(11)
    time_integrate  ...
CubQMCNetG (AbstractStoppingCriterion)
    abs_tol         0.001
    rel_tol         0
    n_init          2^(10)
    n_limit         2^(35)
Keister (AbstractIntegrand)
Gaussian (AbstractTrueMeasure)
    mean            0
    covariance      2^(-1)
    decomp_type     PCA
DigitalNetB2 (AbstractLDDiscreteDistribution)
    d               1
    replications    1
    randomize       LMS DS
    gen_mats_source joe_kuo.6.21201.txt
    order           RADICAL INVERSE
    t               63
    alpha           1
    n_limit         2^(32)
    entropy         7

Vector outputs

>>> f = BoxIntegral(DigitalNetB2(3,seed=7),s=[-1,1])
>>> abs_tol = 1e-3
>>> sc = CubQMCNetG(f,abs_tol=abs_tol,rel_tol=0,check_cone=True)
>>> solution,data = sc.integrate()
>>> solution
array([1.19003352, 0.96068403])
>>> data
Data (Data)
    solution        [1.19  0.961]
    comb_bound_low  [1.189 0.96 ]
    comb_bound_high [1.191 0.962]
    comb_bound_diff [0.001 0.002]
    comb_flags      [ True  True]
    n_total         2^(14)
    n               [16384  1024]
    time_integrate  ...
CubQMCNetG (AbstractStoppingCriterion)
    abs_tol         0.001
    rel_tol         0
    n_init          2^(10)
    n_limit         2^(35)
BoxIntegral (AbstractIntegrand)
    s               [-1  1]
Uniform (AbstractTrueMeasure)
    lower_bound     0
    upper_bound     1
DigitalNetB2 (AbstractLDDiscreteDistribution)
    d               3
    replications    1
    randomize       LMS DS
    gen_mats_source joe_kuo.6.21201.txt
    order           RADICAL INVERSE
    t               63
    alpha           1
    n_limit         2^(32)
    entropy         7
>>> sol3neg1 = -np.pi/4-1/2*np.log(2)+np.log(5+3*np.sqrt(3))
>>> sol31 = np.sqrt(3)/4+1/2*np.log(2+np.sqrt(3))-np.pi/24
>>> true_value = np.array([sol3neg1,sol31])
>>> assert (abs(true_value-solution)<abs_tol).all()

Sensitivity indices 

>>> function = Genz(DigitalNetB2(3,seed=7))
>>> integrand = SensitivityIndices(function)
>>> sc = CubQMCNetG(integrand,abs_tol=5e-4,rel_tol=0,check_cone=True)
>>> solution,data = sc.integrate()
>>> data
Data (Data)
    solution        [[0.02  0.196 0.667]
                     [0.036 0.303 0.782]]
    comb_bound_low  [[0.019 0.195 0.667]
                     [0.035 0.303 0.781]]
    comb_bound_high [[0.02  0.196 0.667]
                     [0.036 0.303 0.782]]
    comb_bound_diff [[0.001 0.    0.001]
                     [0.001 0.001 0.001]]
    comb_flags      [[ True  True  True]
                     [ True  True  True]]
    n_total         2^(16)
    n               [[[16384 65536 65536]
                      [16384 65536 65536]
                      [16384 65536 65536]]

                     [[ 2048 16384 32768]
                      [ 2048 16384 32768]
                      [ 2048 16384 32768]]]
    time_integrate  ...
CubQMCNetG (AbstractStoppingCriterion)
    abs_tol         5.00e-04
    rel_tol         0
    n_init          2^(10)
    n_limit         2^(35)
SensitivityIndices (AbstractIntegrand)
    indices         [[ True False False]
                     [False  True False]
                     [False False  True]]
Uniform (AbstractTrueMeasure)
    lower_bound     0
    upper_bound     1
DigitalNetB2 (AbstractLDDiscreteDistribution)
    d               3
    replications    1
    randomize       LMS DS
    gen_mats_source joe_kuo.6.21201.txt
    order           RADICAL INVERSE
    t               63
    alpha           1
    n_limit         2^(32)
    entropy         7

Control Variates

>>> dnb2 = DigitalNetB2(dimension=4,seed=7)
>>> integrand = CustomFun(
...     true_measure = Uniform(dnb2),
...     g = lambda t: np.stack([
...         1*t[...,0]+2*t[...,0]**2+3*t[...,0]**3,
...         2*t[...,1]+3*t[...,1]**2+4*t[...,1]**3,
...         3*t[...,2]+4*t[...,2]**2+5*t[...,2]**3]),
...     dimension_indv = (3,))
>>> control_variates = [
...     CustomFun(
...         true_measure = Uniform(dnb2),
...         g = lambda t: np.stack([t[...,0],t[...,1],t[...,2]],axis=0),
...         dimension_indv = (3,)),
...     CustomFun(
...         true_measure = Uniform(dnb2),
...         g = lambda t: np.stack([t[...,0]**2,t[...,1]**2,t[...,2]**2],axis=0),
...         dimension_indv = (3,))]
>>> control_variate_means = np.array([[1/2,1/2,1/2],[1/3,1/3,1/3]])
>>> true_value = np.array([23/12,3,49/12])
>>> abs_tol = 1e-6
>>> sc = CubQMCNetG(integrand,abs_tol=abs_tol,rel_tol=0,control_variates=control_variates,control_variate_means=control_variate_means,update_cv_coeffs=False)
>>> solution,data = sc.integrate()
>>> solution
array([1.91666667, 3.        , 4.08333333])
>>> data.n
array([ 8192,  8192, 16384])
>>> assert (np.abs(true_value-solution)<abs_tol).all()
>>> sc = CubQMCNetG(integrand,abs_tol=abs_tol,rel_tol=0,control_variates=control_variates,control_variate_means=control_variate_means,update_cv_coeffs=True)
>>> solution,data = sc.integrate()
>>> solution
array([1.91666667, 3.        , 4.08333333])
>>> data.n
array([16384, 16384, 16384])
>>> assert (np.abs(true_value-solution)<abs_tol).all()

References:

  1. Hickernell, Fred J., and Lluís Antoni Jiménez Rugama.
    "Reliable adaptive cubature using digital sequences."
    Monte Carlo and Quasi-Monte Carlo Methods: MCQMC, Leuven, Belgium, April 2014.
    Springer International Publishing, 2016.

  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.
    http://gailgithub.github.io/GAIL_Dev/.
    https://github.com/GailGithub/GAIL_Dev/blob/master/Algorithms/IntegrationExpectation/cubSobol_g.m.

Parameters:

Name Type Description Default
integrand AbstractIntegrand

The integrand.

 required
abs_tol ndarray

Absolute error tolerance.

 0.01
rel_tol ndarray

Relative error tolerance.

 0.0
n_init int

Initial number of samples.

2 ** 10
n_limit int

Maximum number of samples.

 2 ** 35
error_fun Union[str, callable]

Function mapping the approximate solution, absolute error tolerance, and relative error tolerance to the current error bound.

  • 'EITHER', the default, requires the approximation error must be below either the absolue or relative tolerance. Equivalent to setting 
    error_fun = lambda sv,abs_tol,rel_tol: np.maximum(abs_tol,abs(sv)*rel_tol)
  • 'BOTH' requires the approximation error to be below both the absolue and relative tolerance. Equivalent to setting 
    error_fun = lambda sv,abs_tol,rel_tol: np.minimum(abs_tol,abs(sv)*rel_tol)
'EITHER'
fudge function

Positive function multiplying the finite sum of the Fourier coefficients specified in the cone of functions.

lambda m: 5.0 * 2.0 ** -m
check_cone bool

Whether or not to check if the function falls in the cone.

 False
control_variates list

Integrands to use as control variates, each with the same underlying discrete distribution instance.

 []
control_variate_means ndarray

Means of each control variate.

[]
update_cv_coeffs bool

If set to true, the control variate coefficients are recomputed at each iteration. Otherwise they are estimated once after the initial sampling and then fixed.

 False
Source code in qmcpy/stopping_criterion/cub_qmc_net_g.py 
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
def __init__(self, 
             integrand, 
             abs_tol = 1e-2,
             rel_tol = 0., 
             n_init = 2**10, 
             n_limit = 2**35,
             error_fun = "EITHER",
             fudge = lambda m: 5.*2.**(-m), 
             check_cone = False, 
             control_variates = [], 
             control_variate_means = [], 
             update_cv_coeffs = False,
             ):
    r"""
    Args:
        integrand (AbstractIntegrand): The integrand.
        abs_tol (np.ndarray): Absolute error tolerance.
        rel_tol (np.ndarray): Relative error tolerance.
        n_init (int): Initial number of samples. 
        n_limit (int): Maximum number of samples.
        error_fun (Union[str,callable]): Function mapping the approximate solution, absolute error tolerance, and relative error tolerance to the current error bound.

            - `'EITHER'`, the default, requires the approximation error must be below either the absolue *or* relative tolerance.
                Equivalent to setting
                ```python
                error_fun = lambda sv,abs_tol,rel_tol: np.maximum(abs_tol,abs(sv)*rel_tol)
                ```
            - `'BOTH'` requires the approximation error to be below both the absolue *and* relative tolerance. 
                Equivalent to setting
                ```python
                error_fun = lambda sv,abs_tol,rel_tol: np.minimum(abs_tol,abs(sv)*rel_tol)
                ```
        fudge (function): Positive function multiplying the finite sum of the Fourier coefficients specified in the cone of functions. 
        check_cone (bool): Whether or not to check if the function falls in the cone.
        control_variates (list): Integrands to use as control variates, each with the same underlying discrete distribution instance.
        control_variate_means (np.ndarray): Means of each control variate. 
        update_cv_coeffs (bool): If set to true, the control variate coefficients are recomputed at each iteration. 
            Otherwise they are estimated once after the initial sampling and then fixed.
    """
    super(CubQMCNetG,self).__init__(integrand,abs_tol,rel_tol,n_init,n_limit,fudge,
        check_cone,control_variates,control_variate_means,update_cv_coeffs,
        ptransform = 'none',
        ft = fwht,
        omega = omega_fwht,
        allowed_distribs = [DigitalNetB2],
        cast_complex = False,
        error_fun = error_fun)
    if self.discrete_distrib.order!='RADICAL INVERSE':
        raise ParameterError("CubQMCNet_g requires DigitalNetB2 with 'RADICAL INVERSE' order")

CubQMCLatticeG

Bases: AbstractCubQMCLDG

Quasi-Monte Carlo stopping criterion using rank-1 lattice cubature with guarantees for cones of functions with a predictable decay in the Fourier coefficients.

Examples:

>>> k = Keister(Lattice(seed=7))
>>> sc = CubQMCLatticeG(k,abs_tol=1e-3,rel_tol=0,check_cone=True)
>>> solution,data = sc.integrate()
>>> solution
array(1.38037385)
>>> data
Data (Data)
    solution        1.380
    comb_bound_low  1.380
    comb_bound_high 1.381
    comb_bound_diff 0.001
    comb_flags      1
    n_total         2^(11)
    n               2^(11)
    time_integrate  ...
CubQMCLatticeG (AbstractStoppingCriterion)
    abs_tol         0.001
    rel_tol         0
    n_init          2^(10)
    n_limit         2^(30)
Keister (AbstractIntegrand)
Gaussian (AbstractTrueMeasure)
    mean            0
    covariance      2^(-1)
    decomp_type     PCA
Lattice (AbstractLDDiscreteDistribution)
    d               1
    replications    1
    randomize       SHIFT
    gen_vec_source  kuo.lattice-33002-1024-1048576.9125.txt
    order           RADICAL INVERSE
    n_limit         2^(20)
    entropy         7

Vector outputs

>>> f = BoxIntegral(Lattice(3,seed=11),s=[-1,1])
>>> abs_tol = 1e-3
>>> sc = CubQMCLatticeG(f,abs_tol=abs_tol,rel_tol=0,check_cone=True)
>>> solution,data = sc.integrate()
>>> solution
array([1.18947477, 0.96060862])
>>> data
Data (Data)
    solution        [1.189 0.961]
    comb_bound_low  [1.189 0.96 ]
    comb_bound_high [1.19  0.961]
    comb_bound_diff [0.001 0.001]
    comb_flags      [ True  True]
    n_total         2^(13)
    n               [8192 1024]
    time_integrate  ...
CubQMCLatticeG (AbstractStoppingCriterion)
    abs_tol         0.001
    rel_tol         0
    n_init          2^(10)
    n_limit         2^(30)
BoxIntegral (AbstractIntegrand)
    s               [-1  1]
Uniform (AbstractTrueMeasure)
    lower_bound     0
    upper_bound     1
Lattice (AbstractLDDiscreteDistribution)
    d               3
    replications    1
    randomize       SHIFT
    gen_vec_source  kuo.lattice-33002-1024-1048576.9125.txt
    order           RADICAL INVERSE
    n_limit         2^(20)
    entropy         11
>>> sol3neg1 = -np.pi/4-1/2*np.log(2)+np.log(5+3*np.sqrt(3))
>>> sol31 = np.sqrt(3)/4+1/2*np.log(2+np.sqrt(3))-np.pi/24
>>> true_value = np.array([sol3neg1,sol31])
>>> assert (abs(true_value-solution)<abs_tol).all()

Sensitivity indices 

>>> function = Genz(Lattice(3,seed=7))
>>> integrand = SensitivityIndices(function)
>>> sc = CubQMCLatticeG(integrand,abs_tol=5e-4,rel_tol=0,check_cone=True)
>>> solution,data = sc.integrate()
>>> data
Data (Data)
    solution        [[0.021 0.196 0.667]
                     [0.036 0.303 0.782]]
    comb_bound_low  [[0.02  0.196 0.667]
                     [0.035 0.303 0.781]]
    comb_bound_high [[0.021 0.196 0.667]
                     [0.036 0.303 0.782]]
    comb_bound_diff [[0.001 0.001 0.   ]
                     [0.001 0.001 0.001]]
    comb_flags      [[ True  True  True]
                     [ True  True  True]]
    n_total         2^(16)
    n               [[[16384 32768 65536]
                      [16384 32768 65536]
                      [16384 32768 65536]]

                     [[ 2048 16384 32768]
                      [ 2048 16384 32768]
                      [ 2048 16384 32768]]]
    time_integrate  ...
CubQMCLatticeG (AbstractStoppingCriterion)
    abs_tol         5.00e-04
    rel_tol         0
    n_init          2^(10)
    n_limit         2^(30)
SensitivityIndices (AbstractIntegrand)
    indices         [[ True False False]
                     [False  True False]
                     [False False  True]]
Uniform (AbstractTrueMeasure)
    lower_bound     0
    upper_bound     1
Lattice (AbstractLDDiscreteDistribution)
    d               3
    replications    1
    randomize       SHIFT
    gen_vec_source  kuo.lattice-33002-1024-1048576.9125.txt
    order           RADICAL INVERSE
    n_limit         2^(20)
    entropy         7

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.
    http://gailgithub.github.io/GAIL_Dev/.
    https://github.com/GailGithub/GAIL_Dev/blob/master/Algorithms/IntegrationExpectation/cubLattice_g.m.

Parameters:

Name Type Description Default
integrand AbstractIntegrand

The integrand.

 required
abs_tol ndarray

Absolute error tolerance.

 0.01
rel_tol ndarray

Relative error tolerance.

 0.0
n_init int

Initial number of samples.

2 ** 10
n_limit int

Maximum number of samples.

 2 ** 30
error_fun Union[str, callable]

Function mapping the approximate solution, absolute error tolerance, and relative error tolerance to the current error bound.

  • 'EITHER', the default, requires the approximation error must be below either the absolue or relative tolerance. Equivalent to setting 
    error_fun = lambda sv,abs_tol,rel_tol: np.maximum(abs_tol,abs(sv)*rel_tol)
  • 'BOTH' requires the approximation error to be below both the absolue and relative tolerance. Equivalent to setting 
    error_fun = lambda sv,abs_tol,rel_tol: np.minimum(abs_tol,abs(sv)*rel_tol)
'EITHER'
fudge function

Positive function multiplying the finite sum of the Fourier coefficients specified in the cone of functions.

lambda m: 5.0 * 2.0 ** -m
check_cone bool

Whether or not to check if the function falls in the cone.

 False
ptransform str

Periodization transform, see the options in AbstractIntegrand.f.

 'BAKER'
Source code in qmcpy/stopping_criterion/cub_qmc_lattice_g.py 
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
def __init__(self, 
             integrand, 
             abs_tol = 1e-2, 
             rel_tol = 0., 
             n_init = 2**10,
             n_limit = 2**30,
             error_fun = "EITHER",
             fudge = lambda m: 5.*2.**(-m), 
             check_cone = False,
             ptransform = 'BAKER',
             ):
    r"""
    Args:
        integrand (AbstractIntegrand): The integrand.
        abs_tol (np.ndarray): Absolute error tolerance.
        rel_tol (np.ndarray): Relative error tolerance.
        n_init (int): Initial number of samples. 
        n_limit (int): Maximum number of samples.
        error_fun (Union[str,callable]): Function mapping the approximate solution, absolute error tolerance, and relative error tolerance to the current error bound.

            - `'EITHER'`, the default, requires the approximation error must be below either the absolue *or* relative tolerance.
                Equivalent to setting
                ```python
                error_fun = lambda sv,abs_tol,rel_tol: np.maximum(abs_tol,abs(sv)*rel_tol)
                ```
            - `'BOTH'` requires the approximation error to be below both the absolue *and* relative tolerance. 
                Equivalent to setting
                ```python
                error_fun = lambda sv,abs_tol,rel_tol: np.minimum(abs_tol,abs(sv)*rel_tol)
                ```
        fudge (function): Positive function multiplying the finite sum of the Fourier coefficients specified in the cone of functions. 
        check_cone (bool): Whether or not to check if the function falls in the cone.
        ptransform (str): Periodization transform, see the options in [`AbstractIntegrand.f`][qmcpy.AbstractIntegrand.f].
    """
    super(CubQMCLatticeG,self).__init__(integrand,abs_tol,rel_tol,n_init,n_limit,fudge,check_cone,
        control_variates = [],
        control_variate_means = [],
        update_beta = False,
        ptransform = ptransform,
        ft = fftbr,
        omega = omega_fftbr,
        allowed_distribs = [Lattice],
        cast_complex = True,
        error_fun = error_fun)
    if self.discrete_distrib.order!='RADICAL INVERSE':
        raise ParameterError("CubLattice_g requires Lattice with 'RADICAL INVERSE' order")

CubQMCBayesNetG

Bases: AbstractCubBayesLDG

Quasi-Monte Carlo stopping criterion using fast Bayesian cubature and digital nets with guarantees for Gaussian processes having certain digitally shift invariant kernels.

Examples:

>>> k = Keister(DigitalNetB2(2, seed=123456789))
>>> sc = CubQMCBayesNetG(k,abs_tol=5e-3)
>>> solution,data = sc.integrate()
>>> data
Data (Data)
    solution        1.809
    comb_bound_low  1.807
    comb_bound_high 1.810
    comb_bound_diff 0.003
    comb_flags      1
    n_total         2^(10)
    n               2^(10)
    time_integrate  ...
CubQMCBayesNetG (AbstractStoppingCriterion)
    abs_tol         0.005
    rel_tol         0
    n_init          2^(8)
    n_limit         2^(22)
    order           1
Keister (AbstractIntegrand)
Gaussian (AbstractTrueMeasure)
    mean            0
    covariance      2^(-1)
    decomp_type     PCA
DigitalNetB2 (AbstractLDDiscreteDistribution)
    d               2^(1)
    replications    1
    randomize       LMS DS
    gen_mats_source joe_kuo.6.21201.txt
    order           RADICAL INVERSE
    t               63
    alpha           1
    n_limit         2^(32)
    entropy         123456789

Vector outputs

>>> f = BoxIntegral(DigitalNetB2(3,seed=7),s=[-1,1])
>>> abs_tol = 1e-2
>>> sc = CubQMCBayesNetG(f,abs_tol=abs_tol,rel_tol=0)
>>> solution,data = sc.integrate()
>>> solution
array([1.18750491, 0.96076395])
>>> data
Data (Data)
    solution        [1.188 0.961]
    comb_bound_low  [1.18 0.96]
    comb_bound_high [1.195 0.962]
    comb_bound_diff [0.014 0.002]
    comb_flags      [ True  True]
    n_total         2^(11)
    n               [2048  256]
    time_integrate  ...
CubQMCBayesNetG (AbstractStoppingCriterion)
    abs_tol         0.010
    rel_tol         0
    n_init          2^(8)
    n_limit         2^(22)
    order           1
BoxIntegral (AbstractIntegrand)
    s               [-1  1]
Uniform (AbstractTrueMeasure)
    lower_bound     0
    upper_bound     1
DigitalNetB2 (AbstractLDDiscreteDistribution)
    d               3
    replications    1
    randomize       LMS DS
    gen_mats_source joe_kuo.6.21201.txt
    order           RADICAL INVERSE
    t               63
    alpha           1
    n_limit         2^(32)
    entropy         7
>>> sol3neg1 = -np.pi/4-1/2*np.log(2)+np.log(5+3*np.sqrt(3))
>>> sol31 = np.sqrt(3)/4+1/2*np.log(2+np.sqrt(3))-np.pi/24
>>> true_value = np.array([sol3neg1,sol31])
>>> assert (abs(true_value-solution)<abs_tol).all()

Sensitivity indices 

>>> function = Genz(DigitalNetB2(3,seed=7))
>>> integrand = SensitivityIndices(function)
>>> sc = CubQMCBayesNetG(integrand,abs_tol=5e-2,rel_tol=0)
>>> solution,data = sc.integrate()
>>> data
Data (Data)
    solution        [[0.009 0.194 0.657]
                     [0.036 0.312 0.783]]
    comb_bound_low  [[0.007 0.178 0.636]
                     [0.033 0.297 0.762]]
    comb_bound_high [[0.012 0.209 0.679]
                     [0.038 0.327 0.804]]
    comb_bound_diff [[0.005 0.031 0.043]
                     [0.004 0.03  0.042]]
    comb_flags      [[ True  True  True]
                     [ True  True  True]]
    n_total         2^(8)
    n               [[[256 256 256]
                      [256 256 256]
                      [256 256 256]]

                     [[256 256 256]
                      [256 256 256]
                      [256 256 256]]]
    time_integrate  ...
CubQMCBayesNetG (AbstractStoppingCriterion)
    abs_tol         0.050
    rel_tol         0
    n_init          2^(8)
    n_limit         2^(22)
    order           1
SensitivityIndices (AbstractIntegrand)
    indices         [[ True False False]
                     [False  True False]
                     [False False  True]]
Uniform (AbstractTrueMeasure)
    lower_bound     0
    upper_bound     1
DigitalNetB2 (AbstractLDDiscreteDistribution)
    d               3
    replications    1
    randomize       LMS DS
    gen_mats_source joe_kuo.6.21201.txt
    order           RADICAL INVERSE
    t               63
    alpha           1
    n_limit         2^(32)
    entropy         7

References:

  1. Jagadeeswaran, Rathinavel, and Fred J. Hickernell.
    "Fast automatic Bayesian cubature using Sobol’sampling."
    Advances in Modeling and Simulation: Festschrift for Pierre L'Ecuyer.
    Springer International Publishing, 2022. 301-318.

  2. Jagadeeswaran Rathinavel,
    Fast automatic Bayesian cubature using matching kernels and designs,
    PhD thesis, Illinois Institute of Technology, 2019.

  3. 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.
    http://gailgithub.github.io/GAIL_Dev/.
    https://github.com/GailGithub/GAIL_Dev/blob/master/Algorithms/IntegrationExpectation/cubBayesNet_g.m.

Parameters:

Name Type Description Default
integrand AbstractIntegrand

The integrand.

 required
abs_tol ndarray

Absolute error tolerance.

 0.01
rel_tol ndarray

Relative error tolerance.

 0
n_init int

Initial number of samples.

2 ** 8
n_limit int

Maximum number of samples.

 2 ** 22
error_fun Union[str, callable]

Function mapping the approximate solution, absolute error tolerance, and relative error tolerance to the current error bound.

  • 'EITHER', the default, requires the approximation error must be below either the absolue or relative tolerance. Equivalent to setting 
    error_fun = lambda sv,abs_tol,rel_tol: np.maximum(abs_tol,abs(sv)*rel_tol)
  • 'BOTH' requires the approximation error to be below both the absolue and relative tolerance. Equivalent to setting 
    error_fun = lambda sv,abs_tol,rel_tol: np.minimum(abs_tol,abs(sv)*rel_tol)
'EITHER'
alpha ndarray

Uncertainty level in \((0,1)\).

0.01
errbd_type str

Options are

  • 'MLE': Marginal Log Likelihood.
  • 'GCV': Generalized Cross Validation.
  • 'FULL': Full Bayes.
 'MLE'
Source code in qmcpy/stopping_criterion/cub_qmc_bayes_net_g.py 
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
def __init__(self, 
             integrand, 
             abs_tol = 1e-2, 
             rel_tol = 0,
             n_init = 2**8, 
             n_limit = 2**22, 
             error_fun = "EITHER", 
             alpha = 0.01,
             errbd_type="MLE",
             ):
    r"""
    Args:
        integrand (AbstractIntegrand): The integrand.
        abs_tol (np.ndarray): Absolute error tolerance.
        rel_tol (np.ndarray): Relative error tolerance.
        n_init (int): Initial number of samples. 
        n_limit (int): Maximum number of samples.
        error_fun (Union[str,callable]): Function mapping the approximate solution, absolute error tolerance, and relative error tolerance to the current error bound.

            - `'EITHER'`, the default, requires the approximation error must be below either the absolue *or* relative tolerance.
                Equivalent to setting
                ```python
                error_fun = lambda sv,abs_tol,rel_tol: np.maximum(abs_tol,abs(sv)*rel_tol)
                ```
            - `'BOTH'` requires the approximation error to be below both the absolue *and* relative tolerance. 
                Equivalent to setting
                ```python
                error_fun = lambda sv,abs_tol,rel_tol: np.minimum(abs_tol,abs(sv)*rel_tol)
                ```
        alpha (np.ndarray): Uncertainty level in $(0,1)$. 
        errbd_type (str): Options are 

            - `'MLE'`: Marginal Log Likelihood. 
            - `'GCV'`: Generalized Cross Validation. 
            - `'FULL'`: Full Bayes.
    """
    super(CubQMCBayesNetG, self).__init__(integrand, ft=fwht, omega=omega_fwht,
                                       ptransform=None,
                                       allowed_distribs=[DigitalNetB2],
                                       kernel=self._shift_inv_kernel_digital,
                                       abs_tol=abs_tol, rel_tol=rel_tol,
                                       n_init=n_init, n_limit=n_limit, alpha=alpha, error_fun=error_fun, errbd_type=errbd_type)
    self.order = 1  # Currently supports only order=1
    # private properties
    # Full Bayes - assumes m and s^2 as hyperparameters
    # GCV - Generalized cross validation
    self.kernType = 1  # Type-1:
    self._xfullundtype = np.uint64
    if self.discrete_distrib.order!='RADICAL INVERSE':
        raise ParameterError("CubQMCNet_g requires DigitalNetB2 with 'RADICAL INVERSE' order")

CubQMCBayesLatticeG

Bases: AbstractCubBayesLDG

Quasi-Monte Carlo stopping criterion using fast Bayesian cubature and rank-1 lattices with guarantees for Gaussian processes having certain shift invariant kernels.

Examples:

>>> k = Keister(Lattice(2, seed=123456789))
>>> sc = CubQMCBayesLatticeG(k,abs_tol=1e-4)
>>> solution,data = sc.integrate()
>>> data
Data (Data)
    solution        1.808
    comb_bound_low  1.808
    comb_bound_high 1.808
    comb_bound_diff 3.60e-05
    comb_flags      1
    n_total         2^(10)
    n               2^(10)
    time_integrate  ...
CubQMCBayesLatticeG (AbstractStoppingCriterion)
    abs_tol         1.00e-04
    rel_tol         0
    n_init          2^(8)
    n_limit         2^(22)
    order           2^(1)
Keister (AbstractIntegrand)
Gaussian (AbstractTrueMeasure)
    mean            0
    covariance      2^(-1)
    decomp_type     PCA
Lattice (AbstractLDDiscreteDistribution)
    d               2^(1)
    replications    1
    randomize       SHIFT
    gen_vec_source  kuo.lattice-33002-1024-1048576.9125.txt
    order           RADICAL INVERSE
    n_limit         2^(20)
    entropy         123456789

Vector outputs

>>> f = BoxIntegral(Lattice(3,seed=7),s=[-1,1])
>>> abs_tol = 1e-2
>>> sc = CubQMCBayesLatticeG(f,abs_tol=abs_tol,rel_tol=0)
>>> solution,data = sc.integrate()
>>> solution
array([1.18837601, 0.95984299])
>>> data
Data (Data)
    solution        [1.188 0.96 ]
    comb_bound_low  [1.183 0.95 ]
    comb_bound_high [1.194 0.969]
    comb_bound_diff [0.011 0.019]
    comb_flags      [ True  True]
    n_total         2^(9)
    n               [512 256]
    time_integrate  ...
CubQMCBayesLatticeG (AbstractStoppingCriterion)
    abs_tol         0.010
    rel_tol         0
    n_init          2^(8)
    n_limit         2^(22)
    order           2^(1)
BoxIntegral (AbstractIntegrand)
    s               [-1  1]
Uniform (AbstractTrueMeasure)
    lower_bound     0
    upper_bound     1
Lattice (AbstractLDDiscreteDistribution)
    d               3
    replications    1
    randomize       SHIFT
    gen_vec_source  kuo.lattice-33002-1024-1048576.9125.txt
    order           RADICAL INVERSE
    n_limit         2^(20)
    entropy         7
>>> sol3neg1 = -np.pi/4-1/2*np.log(2)+np.log(5+3*np.sqrt(3))
>>> sol31 = np.sqrt(3)/4+1/2*np.log(2+np.sqrt(3))-np.pi/24
>>> true_value = np.array([sol3neg1,sol31])
>>> assert (abs(true_value-solution)<abs_tol).all()

Sensitivity indices 

>>> function = Genz(Lattice(3,seed=7))
>>> integrand = SensitivityIndices(function)
>>> sc = CubQMCBayesLatticeG(integrand,abs_tol=5e-2,rel_tol=0)
>>> solution,data = sc.integrate()
>>> data
Data (Data)
    solution        [[0.057 0.131 0.269]
                     [0.386 0.523 0.741]]
    comb_bound_low  [[0.048 0.12  0.256]
                     [0.377 0.511 0.721]]
    comb_bound_high [[0.066 0.141 0.283]
                     [0.395 0.535 0.762]]
    comb_bound_diff [[0.018 0.021 0.027]
                     [0.018 0.024 0.04 ]]
    comb_flags      [[ True  True  True]
                     [ True  True  True]]
    n_total         2^(10)
    n               [[[1024 1024 1024]
                      [1024 1024 1024]
                      [1024 1024 1024]]

                     [[1024 1024 1024]
                      [1024 1024 1024]
                      [1024 1024 1024]]]
    time_integrate  ...
CubQMCBayesLatticeG (AbstractStoppingCriterion)
    abs_tol         0.050
    rel_tol         0
    n_init          2^(8)
    n_limit         2^(22)
    order           2^(1)
SensitivityIndices (AbstractIntegrand)
    indices         [[ True False False]
                     [False  True False]
                     [False False  True]]
Uniform (AbstractTrueMeasure)
    lower_bound     0
    upper_bound     1
Lattice (AbstractLDDiscreteDistribution)
    d               3
    replications    1
    randomize       SHIFT
    gen_vec_source  kuo.lattice-33002-1024-1048576.9125.txt
    order           RADICAL INVERSE
    n_limit         2^(20)
    entropy         7

References:

  1. Jagadeeswaran, Rathinavel, and Fred J. Hickernell.
    "Fast automatic Bayesian cubature using lattice sampling."
    Statistics and Computing 29.6 (2019): 1215-1229.

  2. Jagadeeswaran Rathinavel and Fred J. Hickernell,
    Fast automatic Bayesian cubature using lattice sampling.
    Stat Comput 29, 1215-1229 (2019).
    Available from Springer https://doi.org/10.1007/s11222-019-09895-9.

  3. 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.
    http://gailgithub.github.io/GAIL_Dev/.
    https://github.com/GailGithub/GAIL_Dev/blob/master/Algorithms/IntegrationExpectation/cubBayesLattice_g.m.

Parameters:

Name Type Description Default
integrand AbstractIntegrand

The integrand.

 required
abs_tol ndarray

Absolute error tolerance.

 0.01
rel_tol ndarray

Relative error tolerance.

 0
n_init int

Initial number of samples.

2 ** 8
n_limit int

Maximum number of samples.

 2 ** 22
error_fun Union[str, callable]

Function mapping the approximate solution, absolute error tolerance, and relative error tolerance to the current error bound.

  • 'EITHER', the default, requires the approximation error must be below either the absolue or relative tolerance. Equivalent to setting 
    error_fun = lambda sv,abs_tol,rel_tol: np.maximum(abs_tol,abs(sv)*rel_tol)
  • 'BOTH' requires the approximation error to be below both the absolue and relative tolerance. Equivalent to setting 
    error_fun = lambda sv,abs_tol,rel_tol: np.minimum(abs_tol,abs(sv)*rel_tol)
'EITHER'
alpha ndarray

Uncertainty level in \((0,1)\).

0.01
ptransform str

Periodization transform, see the options in AbstractIntegrand.f.

 'C1SIN'
errbd_type str

Options are

  • 'MLE': Marginal Log Likelihood.
  • 'GCV': Generalized Cross Validation.
  • 'FULL': Full Bayes.
 'MLE'
order int

Bernoulli kernel's order. If zero, choose order automatically

 2
Source code in qmcpy/stopping_criterion/cub_qmc_bayes_lattice_g.py 
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
def __init__(self, 
             integrand, 
             abs_tol = 1e-2, 
             rel_tol = 0,
             n_init = 2**8, 
             n_limit = 2**22, 
             error_fun="EITHER", 
             alpha = 0.01, 
             ptransform = 'C1SIN',
             errbd_type = "MLE", 
             order = 2,
             ):
    r"""
    Args:
        integrand (AbstractIntegrand): The integrand.
        abs_tol (np.ndarray): Absolute error tolerance.
        rel_tol (np.ndarray): Relative error tolerance.
        n_init (int): Initial number of samples. 
        n_limit (int): Maximum number of samples.
        error_fun (Union[str,callable]): Function mapping the approximate solution, absolute error tolerance, and relative error tolerance to the current error bound.

            - `'EITHER'`, the default, requires the approximation error must be below either the absolue *or* relative tolerance.
                Equivalent to setting
                ```python
                error_fun = lambda sv,abs_tol,rel_tol: np.maximum(abs_tol,abs(sv)*rel_tol)
                ```
            - `'BOTH'` requires the approximation error to be below both the absolue *and* relative tolerance. 
                Equivalent to setting
                ```python
                error_fun = lambda sv,abs_tol,rel_tol: np.minimum(abs_tol,abs(sv)*rel_tol)
                ```
        alpha (np.ndarray): Uncertainty level in $(0,1)$. 
        ptransform (str): Periodization transform, see the options in [`AbstractIntegrand.f`][qmcpy.AbstractIntegrand.f].
        errbd_type (str): Options are 

            - `'MLE'`: Marginal Log Likelihood. 
            - `'GCV'`: Generalized Cross Validation. 
            - `'FULL'`: Full Bayes.
        order (int): Bernoulli kernel's order. If zero, choose order automatically
    """
    super(CubQMCBayesLatticeG, self).__init__(integrand, ft=fftbr, omega=omega_fftbr,
                                           ptransform=ptransform,
                                           allowed_distribs=[Lattice],
                                           kernel=self._shift_inv_kernel,
                                           abs_tol=abs_tol, rel_tol=rel_tol,
                                           n_init=n_init, n_limit=n_limit, alpha=alpha, error_fun=error_fun, errbd_type=errbd_type)
    self.order = order  # Bernoulli kernel's order. If zero, choose order automatically
    # private properties
    # full_Bayes - Full Bayes - assumes m and s^2 as hyperparameters
    # GCV - Generalized cross validation
    self.kernType = 1  # Type-1: Bernoulli polynomial based algebraic convergence, Type-2: Truncated series
    self._xfullundtype = float
    if self.discrete_distrib.order!='RADICAL INVERSE':
        raise ParameterError("CubLattice_g requires Lattice with 'RADICAL INVERSE' order")

CubQMCRepStudentT

Bases: AbstractStoppingCriterion

Quasi-Monte Carlo stopping criterion based on Student's \(t\)-distribution for multiple replications.

Examples:

>>> k = Keister(DigitalNetB2(seed=7,replications=25))
>>> sc = CubQMCRepStudentT(k,abs_tol=1e-3,rel_tol=0)
>>> solution,data = sc.integrate()
>>> solution
array(1.3803927)
>>> data
Data (Data)
    solution        1.380
    comb_bound_low  1.380
    comb_bound_high 1.381
    comb_bound_diff 0.002
    comb_flags      1
    n_total         6400
    n               6400
    n_rep           2^(8)
    time_integrate  ...
CubQMCRepStudentT (AbstractStoppingCriterion)
    inflate         1
    alpha           0.010
    abs_tol         0.001
    rel_tol         0
    n_init          2^(8)
    n_limit         2^(30)
Keister (AbstractIntegrand)
Gaussian (AbstractTrueMeasure)
    mean            0
    covariance      2^(-1)
    decomp_type     PCA
DigitalNetB2 (AbstractLDDiscreteDistribution)
    d               1
    replications    25
    randomize       LMS DS
    gen_mats_source joe_kuo.6.21201.txt
    order           RADICAL INVERSE
    t               63
    alpha           1
    n_limit         2^(32)
    entropy         7

Vector outputs

>>> f = BoxIntegral(DigitalNetB2(3,seed=7,replications=25),s=[-1,1])
>>> abs_tol = 1e-3
>>> sc = CubQMCRepStudentT(f,abs_tol=abs_tol,rel_tol=0)
>>> solution,data = sc.integrate()
>>> solution
array([1.19025707, 0.96062762])
>>> data
Data (Data)
    solution        [1.19  0.961]
    comb_bound_low  [1.19  0.961]
    comb_bound_high [1.191 0.961]
    comb_bound_diff [0.001 0.   ]
    comb_flags      [ True  True]
    n_total         409600
    n               [409600   6400]
    n_rep           [16384   256]
    time_integrate  ...
CubQMCRepStudentT (AbstractStoppingCriterion)
    inflate         1
    alpha           0.010
    abs_tol         0.001
    rel_tol         0
    n_init          2^(8)
    n_limit         2^(30)
BoxIntegral (AbstractIntegrand)
    s               [-1  1]
Uniform (AbstractTrueMeasure)
    lower_bound     0
    upper_bound     1
DigitalNetB2 (AbstractLDDiscreteDistribution)
    d               3
    replications    25
    randomize       LMS DS
    gen_mats_source joe_kuo.6.21201.txt
    order           RADICAL INVERSE
    t               63
    alpha           1
    n_limit         2^(32)
    entropy         7
>>> sol3neg1 = -np.pi/4-1/2*np.log(2)+np.log(5+3*np.sqrt(3))
>>> sol31 = np.sqrt(3)/4+1/2*np.log(2+np.sqrt(3))-np.pi/24
>>> true_value = np.array([sol3neg1,sol31])
>>> assert (abs(true_value-solution)<abs_tol).all()

Sensitivity indices 

>>> function = Genz(DigitalNetB2(3,seed=7,replications=25))
>>> integrand = SensitivityIndices(function)
>>> sc = CubQMCRepStudentT(integrand,abs_tol=5e-4,rel_tol=0)
>>> solution,data = sc.integrate()
>>> data
Data (Data)
    solution        [[0.02  0.196 0.667]
                     [0.036 0.303 0.782]]
    comb_bound_low  [[0.019 0.195 0.667]
                     [0.035 0.303 0.781]]
    comb_bound_high [[0.02  0.196 0.667]
                     [0.036 0.303 0.782]]
    comb_bound_diff [[0.001 0.001 0.001]
                     [0.001 0.001 0.001]]
    comb_flags      [[ True  True  True]
                     [ True  True  True]]
    n_total         204800
    n               [[[102400 204800 204800]
                      [102400 204800 204800]
                      [102400 204800 204800]]

                     [[ 12800 102400 102400]
                      [ 12800 102400 102400]
                      [ 12800 102400 102400]]]
    n_rep           [[[4096 8192 8192]
                      [4096 8192 8192]
                      [4096 8192 8192]]

                     [[ 512 4096 4096]
                      [ 512 4096 4096]
                      [ 512 4096 4096]]]
    time_integrate  ...
CubQMCRepStudentT (AbstractStoppingCriterion)
    inflate         1
    alpha           0.010
    abs_tol         5.00e-04
    rel_tol         0
    n_init          2^(8)
    n_limit         2^(30)
SensitivityIndices (AbstractIntegrand)
    indices         [[ True False False]
                     [False  True False]
                     [False False  True]]
Uniform (AbstractTrueMeasure)
    lower_bound     0
    upper_bound     1
DigitalNetB2 (AbstractLDDiscreteDistribution)
    d               3
    replications    25
    randomize       LMS DS
    gen_mats_source joe_kuo.6.21201.txt
    order           RADICAL INVERSE
    t               63
    alpha           1
    n_limit         2^(32)
    entropy         7

References:

  1. Art B. Owen. "Practical Quasi-Monte Carlo Integration." 2023.
    https://artowen.su.domains/mc/.

  2. Pierre l’Ecuyer et al.
    "Confidence intervals for randomized quasi-Monte Carlo estimators."
    2023 Winter Simulation Conference (WSC). IEEE, 2023.
    https://ieeexplore.ieee.org/stamp/stamp.jsp?arnumber=10408613.

Parameters:

Name Type Description Default
integrand AbstractIntegrand

The integrand.

 required
abs_tol ndarray

Absolute error tolerance.

 0.01
rel_tol ndarray

Relative error tolerance.

 0.0
n_init int

Initial number of samples.

256.0
n_limit int

Maximum number of samples.

 2 ** 30
error_fun Union[str, callable]

Function mapping the approximate solution, absolute error tolerance, and relative error tolerance to the current error bound.

  • 'EITHER', the default, requires the approximation error must be below either the absolue or relative tolerance. Equivalent to setting 
    error_fun = lambda sv,abs_tol,rel_tol: np.maximum(abs_tol,abs(sv)*rel_tol)
  • 'BOTH' requires the approximation error to be below both the absolue and relative tolerance. Equivalent to setting 
    error_fun = lambda sv,abs_tol,rel_tol: np.minimum(abs_tol,abs(sv)*rel_tol)
'EITHER'
inflate float

Inflation factor \(\geq 1\) to multiply by the variance estimate to make it more conservative.

 1
alpha ndarray

Uncertainty level in \((0,1)\).

 0.01
Source code in qmcpy/stopping_criterion/cub_qmc_rep_student_t.py 
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
def __init__(self,
             integrand, 
             abs_tol = 1e-2,
             rel_tol = 0.,
             n_init = 256.,
             n_limit = 2**30,
             error_fun = "EITHER",
             inflate = 1,
             alpha = 0.01, 
             ):
    r"""
    Args:
        integrand (AbstractIntegrand): The integrand.
        abs_tol (np.ndarray): Absolute error tolerance.
        rel_tol (np.ndarray): Relative error tolerance.
        n_init (int): Initial number of samples. 
        n_limit (int): Maximum number of samples.
        error_fun (Union[str,callable]): Function mapping the approximate solution, absolute error tolerance, and relative error tolerance to the current error bound.

            - `'EITHER'`, the default, requires the approximation error must be below either the absolue *or* relative tolerance.
                Equivalent to setting
                ```python
                error_fun = lambda sv,abs_tol,rel_tol: np.maximum(abs_tol,abs(sv)*rel_tol)
                ```
            - `'BOTH'` requires the approximation error to be below both the absolue *and* relative tolerance. 
                Equivalent to setting
                ```python
                error_fun = lambda sv,abs_tol,rel_tol: np.minimum(abs_tol,abs(sv)*rel_tol)
                ```
        inflate (float): Inflation factor $\geq 1$ to multiply by the variance estimate to make it more conservative.
        alpha (np.ndarray): Uncertainty level in $(0,1)$. 
    """
    self.parameters = ['inflate','alpha','abs_tol','rel_tol','n_init','n_limit']
    # Input Checks
    if np.log2(n_init)%1!=0:
        warnings.warn('n_init must be a power of two. Using n_init = 2**5',ParameterWarning)
        n_init = 2**5
    if np.log2(n_limit)%1!=0:
        warnings.warn('n_init must be a power of two. Using n_limit = 2**30',ParameterWarning)
        n_limit = 2**30
    # Set Attributes
    self.n_init = int(n_init)
    self.n_limit = int(n_limit)
    assert isinstance(error_fun,str) or callable(error_fun)
    if isinstance(error_fun,str):
        if error_fun.upper()=="EITHER":
            error_fun = lambda sv,abs_tol,rel_tol: np.maximum(abs_tol,abs(sv)*rel_tol)
        elif error_fun.upper()=="BOTH":
            error_fun = lambda sv,abs_tol,rel_tol: np.minimum(abs_tol,abs(sv)*rel_tol)
        else:
            raise ParameterError("str error_fun must be 'EITHER' or 'BOTH'")
    self.error_fun = error_fun
    self.alpha = alpha
    self.inflate = float(inflate)
    assert self.inflate>=1
    assert 0<self.alpha<1
    # QMCPy Objs
    self.integrand = integrand
    self.true_measure = self.integrand.true_measure
    self.discrete_distrib = self.integrand.discrete_distrib
    super(CubQMCRepStudentT,self).__init__(allowed_distribs=[AbstractLDDiscreteDistribution],allow_vectorized_integrals=True)
    assert self.integrand.discrete_distrib.replications>1, "Require the discrete distribution has replications>1"
    assert self.integrand.discrete_distrib.randomize!="FALSE", "Require discrete distribution is randomized"
    self.alphas_indv,identity_dependency = self._compute_indv_alphas(np.full(self.integrand.d_comb,self.alpha))
    self.set_tolerance(abs_tol,rel_tol)
    self.t_star = -t.ppf(self.alphas_indv/2,df=self.integrand.discrete_distrib.replications-1)

IID MC

CubMCCLTVec

Bases: AbstractStoppingCriterion

IID Monte Carlo stopping criterion stopping criterion based on the Central Limit Theorem with doubling sample sizes.

Examples:

>>> k = Keister(IIDStdUniform(seed=7))
>>> sc = CubMCCLTVec(k,abs_tol=5e-2,rel_tol=0)
>>> solution,data = sc.integrate()
>>> solution
array(1.38366791)
>>> data
Data (Data)
    solution        1.384
    comb_bound_low  1.343
    comb_bound_high 1.424
    comb_bound_diff 0.080
    comb_flags      1
    n_total         1024
    n               2^(10)
    time_integrate  ...
CubMCCLTVec (AbstractStoppingCriterion)
    inflate         1
    alpha           0.010
    abs_tol         0.050
    rel_tol         0
    n_init          2^(8)
    n_limit         2^(30)
Keister (AbstractIntegrand)
Gaussian (AbstractTrueMeasure)
    mean            0
    covariance      2^(-1)
    decomp_type     PCA
IIDStdUniform (AbstractIIDDiscreteDistribution)
    d               1
    replications    1
    entropy         7

Vector outputs

>>> f = BoxIntegral(IIDStdUniform(3,seed=7),s=[-1,1])
>>> abs_tol = 2.5e-2
>>> sc = CubMCCLTVec(f,abs_tol=abs_tol,rel_tol=0)
>>> solution,data = sc.integrate()
>>> solution
array([1.18448043, 0.95435347])
>>> data
Data (Data)
    solution        [1.184 0.954]
    comb_bound_low  [1.165 0.932]
    comb_bound_high [1.203 0.977]
    comb_bound_diff [0.038 0.045]
    comb_flags      [ True  True]
    n_total         8192
    n               [8192 1024]
    time_integrate  ...
CubMCCLTVec (AbstractStoppingCriterion)
    inflate         1
    alpha           0.010
    abs_tol         0.025
    rel_tol         0
    n_init          2^(8)
    n_limit         2^(30)
BoxIntegral (AbstractIntegrand)
    s               [-1  1]
Uniform (AbstractTrueMeasure)
    lower_bound     0
    upper_bound     1
IIDStdUniform (AbstractIIDDiscreteDistribution)
    d               3
    replications    1
    entropy         7
>>> sol3neg1 = -np.pi/4-1/2*np.log(2)+np.log(5+3*np.sqrt(3))
>>> sol31 = np.sqrt(3)/4+1/2*np.log(2+np.sqrt(3))-np.pi/24
>>> true_value = np.array([sol3neg1,sol31])
>>> assert (abs(true_value-solution)<abs_tol).all()

Sensitivity indices 

>>> function = Genz(IIDStdUniform(3,seed=7))
>>> integrand = SensitivityIndices(function)
>>> sc = CubMCCLTVec(integrand,abs_tol=2.5e-2,rel_tol=0)
>>> solution,data = sc.integrate()
>>> data
Data (Data)
    solution        [[0.024 0.203 0.662]
                     [0.044 0.308 0.78 ]]
    comb_bound_low  [[0.006 0.186 0.644]
                     [0.02  0.286 0.761]]
    comb_bound_high [[0.042 0.221 0.681]
                     [0.067 0.329 0.798]]
    comb_bound_diff [[0.036 0.035 0.037]
                     [0.047 0.043 0.037]]
    comb_flags      [[ True  True  True]
                     [ True  True  True]]
    n_total         262144
    n               [[[  4096  65536 262144]
                      [  4096  65536 262144]
                      [  4096  65536 262144]]

                     [[   512  32768 262144]
                      [   512  32768 262144]
                      [   512  32768 262144]]]
    time_integrate  ...
CubMCCLTVec (AbstractStoppingCriterion)
    inflate         1
    alpha           0.010
    abs_tol         0.025
    rel_tol         0
    n_init          2^(8)
    n_limit         2^(30)
SensitivityIndices (AbstractIntegrand)
    indices         [[ True False False]
                     [False  True False]
                     [False False  True]]
Uniform (AbstractTrueMeasure)
    lower_bound     0
    upper_bound     1
IIDStdUniform (AbstractIIDDiscreteDistribution)
    d               3
    replications    1
    entropy         7

Parameters:

Name Type Description Default
integrand AbstractIntegrand

The integrand.

 required
abs_tol ndarray

Absolute error tolerance.

 0.01
rel_tol ndarray

Relative error tolerance.

 0.0
n_init int

Initial number of samples.

256.0
n_limit int

Maximum number of samples.

 2 ** 30
error_fun Union[str, callable]

Function mapping the approximate solution, absolute error tolerance, and relative error tolerance to the current error bound.

  • 'EITHER', the default, requires the approximation error must be below either the absolue or relative tolerance. Equivalent to setting 
    error_fun = lambda sv,abs_tol,rel_tol: np.maximum(abs_tol,abs(sv)*rel_tol)
  • 'BOTH' requires the approximation error to be below both the absolue and relative tolerance. Equivalent to setting 
    error_fun = lambda sv,abs_tol,rel_tol: np.minimum(abs_tol,abs(sv)*rel_tol)
'EITHER'
inflate float

Inflation factor \(\geq 1\) to multiply by the variance estimate to make it more conservative.

 1
alpha ndarray

Uncertainty level in \((0,1)\).

 0.01
Source code in qmcpy/stopping_criterion/cub_mc_clt_vec.py 
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
def __init__(self,
             integrand, 
             abs_tol = 1e-2,
             rel_tol = 0.,
             n_init = 256.,
             n_limit = 2**30,
             error_fun = "EITHER",
             inflate = 1,
             alpha = 0.01,
             ):
    r"""
    Args:
        integrand (AbstractIntegrand): The integrand.
        abs_tol (np.ndarray): Absolute error tolerance.
        rel_tol (np.ndarray): Relative error tolerance.
        n_init (int): Initial number of samples. 
        n_limit (int): Maximum number of samples.
        error_fun (Union[str,callable]): Function mapping the approximate solution, absolute error tolerance, and relative error tolerance to the current error bound.

            - `'EITHER'`, the default, requires the approximation error must be below either the absolue *or* relative tolerance.
                Equivalent to setting
                ```python
                error_fun = lambda sv,abs_tol,rel_tol: np.maximum(abs_tol,abs(sv)*rel_tol)
                ```
            - `'BOTH'` requires the approximation error to be below both the absolue *and* relative tolerance. 
                Equivalent to setting
                ```python
                error_fun = lambda sv,abs_tol,rel_tol: np.minimum(abs_tol,abs(sv)*rel_tol)
                ```
        inflate (float): Inflation factor $\geq 1$ to multiply by the variance estimate to make it more conservative.
        alpha (np.ndarray): Uncertainty level in $(0,1)$. 
    """
    self.parameters = ['inflate','alpha','abs_tol','rel_tol','n_init','n_limit']
    # Input Checks
    if np.log2(n_init)%1!=0:
        warnings.warn('n_init must be a power of two. Using n_init = 2**5',ParameterWarning)
        n_init = 2**5
    if np.log2(n_limit)%1!=0:
        warnings.warn('n_init must be a power of two. Using n_limit = 2**30',ParameterWarning)
        n_limit = 2**30
    # Set Attributes
    self.n_init = int(n_init)
    self.n_limit = int(n_limit)
    assert isinstance(error_fun,str) or callable(error_fun)
    if isinstance(error_fun,str):
        if error_fun.upper()=="EITHER":
            error_fun = lambda sv,abs_tol,rel_tol: np.maximum(abs_tol,abs(sv)*rel_tol)
        elif error_fun.upper()=="BOTH":
            error_fun = lambda sv,abs_tol,rel_tol: np.minimum(abs_tol,abs(sv)*rel_tol)
        else:
            raise ParameterError("str error_fun must be 'EITHER' or 'BOTH'")
    self.error_fun = error_fun
    self.alpha = alpha
    self.inflate = float(inflate)
    assert self.inflate>=1
    # QMCPy Objs
    self.integrand = integrand
    self.true_measure = self.integrand.true_measure
    self.discrete_distrib = self.integrand.discrete_distrib
    super(CubMCCLTVec,self).__init__(allowed_distribs=[AbstractIIDDiscreteDistribution],allow_vectorized_integrals=True)
    assert self.integrand.discrete_distrib.no_replications==True, "Require the discrete distribution has replications=None"
    self.alphas_indv,identity_dependency = self._compute_indv_alphas(np.full(self.integrand.d_comb,self.alpha))
    self.set_tolerance(abs_tol,rel_tol)
    self.z_star = -norm.ppf(self.alphas_indv/2)

CubMCG

Bases: AbstractStoppingCriterion

IID Monte Carlo stopping criterion using Berry-Esseen inequalities in a two step method with guarantees for functions with bounded kurtosis.

Examples:

>>> ao = FinancialOption(IIDStdUniform(52,seed=7))
>>> sc = CubMCG(ao,abs_tol=.05)
>>> solution,data = sc.integrate()
>>> data
Data (Data)
    solution        1.779
    bound_low       1.729
    bound_high      1.829
    bound_diff      0.100
    n_total         112314
    time_integrate  ...
CubMCG (AbstractStoppingCriterion)
    abs_tol         0.050
    rel_tol         0
    n_init          2^(10)
    n_limit         2^(30)
    inflate         1.200
    alpha           0.010
    kurtmax         1.478
FinancialOption (AbstractIntegrand)
    option          ASIAN
    call_put        CALL
    volatility      2^(-1)
    start_price     30
    strike_price    35
    interest_rate   0
    t_final         1
    asian_mean      ARITHMETIC
BrownianMotion (AbstractTrueMeasure)
    time_vec        [0.019 0.038 0.058 ... 0.962 0.981 1.   ]
    drift           0
    mean            [0. 0. 0. ... 0. 0. 0.]
    covariance      [[0.019 0.019 0.019 ... 0.019 0.019 0.019]
                     [0.019 0.038 0.038 ... 0.038 0.038 0.038]
                     [0.019 0.038 0.058 ... 0.058 0.058 0.058]
                     ...
                     [0.019 0.038 0.058 ... 0.962 0.962 0.962]
                     [0.019 0.038 0.058 ... 0.962 0.981 0.981]
                     [0.019 0.038 0.058 ... 0.962 0.981 1.   ]]
    decomp_type     PCA
IIDStdUniform (AbstractIIDDiscreteDistribution)
    d               52
    replications    1
    entropy         7

Control variates 

>>> iid = IIDStdUniform(52,seed=7)
>>> ao = FinancialOption(iid,option="ASIAN")
>>> eo = FinancialOption(iid,option="EUROPEAN")
>>> lin0 = Linear0(iid)
>>> sc = CubMCG(ao,abs_tol=.05,
...     control_variates = [eo,lin0],
...     control_variate_means = [eo.get_exact_value(),0])
>>> solution,data = sc.integrate()
>>> data
Data (Data)
    solution        1.787
    bound_low       1.737
    bound_high      1.837
    bound_diff      0.100
    n_total         52147
    time_integrate  ...
CubMCG (AbstractStoppingCriterion)
    abs_tol         0.050
    rel_tol         0
    n_init          2^(10)
    n_limit         2^(30)
    inflate         1.200
    alpha           0.010
    kurtmax         1.478
    cv              [FinancialOption (AbstractIntegrand)
                         option          EUROPEAN
                         call_put        CALL
                         volatility      2^(-1)
                         start_price     30
                         strike_price    35
                         interest_rate   0
                         t_final         1               Linear0 (AbstractIntegrand)]
    cv_mu           [4.211 0.   ]
FinancialOption (AbstractIntegrand)
    option          ASIAN
    call_put        CALL
    volatility      2^(-1)
    start_price     30
    strike_price    35
    interest_rate   0
    t_final         1
    asian_mean      ARITHMETIC
BrownianMotion (AbstractTrueMeasure)
    time_vec        [0.019 0.038 0.058 ... 0.962 0.981 1.   ]
    drift           0
    mean            [0. 0. 0. ... 0. 0. 0.]
    covariance      [[0.019 0.019 0.019 ... 0.019 0.019 0.019]
                     [0.019 0.038 0.038 ... 0.038 0.038 0.038]
                     [0.019 0.038 0.058 ... 0.058 0.058 0.058]
                     ...
                     [0.019 0.038 0.058 ... 0.962 0.962 0.962]
                     [0.019 0.038 0.058 ... 0.962 0.981 0.981]
                     [0.019 0.038 0.058 ... 0.962 0.981 1.   ]]
    decomp_type     PCA
IIDStdUniform (AbstractIIDDiscreteDistribution)
    d               52
    replications    1
    entropy         7

Relative tolerance 

>>> ao = FinancialOption(IIDStdUniform(52,seed=7))
>>> sc = CubMCG(ao,abs_tol=1e-3,rel_tol=5e-2)
>>> solution,data = sc.integrate()
>>> data
Data (Data)
    solution        1.743
    bound_low       1.661
    bound_high      1.825
    bound_diff      0.164
    n_total         27503
    time_integrate  ...
CubMCG (AbstractStoppingCriterion)
    abs_tol         0.001
    rel_tol         0.050
    n_init          2^(10)
    n_limit         2^(30)
    inflate         1.200
    alpha           0.010
    kurtmax         1.478
FinancialOption (AbstractIntegrand)
    option          ASIAN
    call_put        CALL
    volatility      2^(-1)
    start_price     30
    strike_price    35
    interest_rate   0
    t_final         1
    asian_mean      ARITHMETIC
BrownianMotion (AbstractTrueMeasure)
    time_vec        [0.019 0.038 0.058 ... 0.962 0.981 1.   ]
    drift           0
    mean            [0. 0. 0. ... 0. 0. 0.]
    covariance      [[0.019 0.019 0.019 ... 0.019 0.019 0.019]
                     [0.019 0.038 0.038 ... 0.038 0.038 0.038]
                     [0.019 0.038 0.058 ... 0.058 0.058 0.058]
                     ...
                     [0.019 0.038 0.058 ... 0.962 0.962 0.962]
                     [0.019 0.038 0.058 ... 0.962 0.981 0.981]
                     [0.019 0.038 0.058 ... 0.962 0.981 1.   ]]
    decomp_type     PCA
IIDStdUniform (AbstractIIDDiscreteDistribution)
    d               52
    replications    1
    entropy         7

Relative tolerance and control variates 

>>> iid = IIDStdUniform(52,seed=7)
>>> ao = FinancialOption(iid,option="ASIAN")
>>> eo = FinancialOption(iid,option="EUROPEAN")
>>> lin0 = Linear0(iid)
>>> sc = CubMCG(ao,abs_tol=1e-3,rel_tol=5e-2,
...     control_variates = [eo,lin0],
...     control_variate_means = [eo.get_exact_value(),0])
>>> solution,data = sc.integrate()
>>> data
Data (Data)
    solution        1.776
    bound_low       1.692
    bound_high      1.859
    bound_diff      0.167
    n_total         12074
    time_integrate  ...
CubMCG (AbstractStoppingCriterion)
    abs_tol         0.001
    rel_tol         0.050
    n_init          2^(10)
    n_limit         2^(30)
    inflate         1.200
    alpha           0.010
    kurtmax         1.478
    cv              [FinancialOption (AbstractIntegrand)
                         option          EUROPEAN
                         call_put        CALL
                         volatility      2^(-1)
                         start_price     30
                         strike_price    35
                         interest_rate   0
                         t_final         1               Linear0 (AbstractIntegrand)]
    cv_mu           [4.211 0.   ]
FinancialOption (AbstractIntegrand)
    option          ASIAN
    call_put        CALL
    volatility      2^(-1)
    start_price     30
    strike_price    35
    interest_rate   0
    t_final         1
    asian_mean      ARITHMETIC
BrownianMotion (AbstractTrueMeasure)
    time_vec        [0.019 0.038 0.058 ... 0.962 0.981 1.   ]
    drift           0
    mean            [0. 0. 0. ... 0. 0. 0.]
    covariance      [[0.019 0.019 0.019 ... 0.019 0.019 0.019]
                     [0.019 0.038 0.038 ... 0.038 0.038 0.038]
                     [0.019 0.038 0.058 ... 0.058 0.058 0.058]
                     ...
                     [0.019 0.038 0.058 ... 0.962 0.962 0.962]
                     [0.019 0.038 0.058 ... 0.962 0.981 0.981]
                     [0.019 0.038 0.058 ... 0.962 0.981 1.   ]]
    decomp_type     PCA
IIDStdUniform (AbstractIIDDiscreteDistribution)
    d               52
    replications    1
    entropy         7

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.
    http://gailgithub.github.io/GAIL_Dev/.
    https://github.com/GailGithub/GAIL_Dev/blob/master/Algorithms/IntegrationExpectation/meanMC_g.m.

Parameters:

Name Type Description Default
integrand AbstractIntegrand

The integrand.

 required
abs_tol ndarray

Absolute error tolerance.

 0.01
rel_tol ndarray

Relative error tolerance.

 0.0
n_init int

Initial number of samples.

1024
n_limit int

Maximum number of samples.

 2 ** 30
inflate float

Inflation factor \(\geq 1\) to multiply by the variance estimate to make it more conservative.

 1.2
alpha ndarray

Uncertainty level in \((0,1)\).

0.01
control_variates list

Integrands to use as control variates, each with the same underlying discrete distribution instance.

 []
control_variate_means ndarray

Means of each control variate.

 []
Source code in qmcpy/stopping_criterion/cub_mc_g.py 
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
def __init__(self, 
             integrand, 
             abs_tol = 1e-2, 
             rel_tol = 0., 
             n_init = 1024, 
             n_limit = 2**30,
             inflate = 1.2, 
             alpha = 0.01, 
             control_variates = [], 
             control_variate_means = [],
             ):
    r"""
    Args:
        integrand (AbstractIntegrand): The integrand.
        abs_tol (np.ndarray): Absolute error tolerance.
        rel_tol (np.ndarray): Relative error tolerance.
        n_init (int): Initial number of samples. 
        n_limit (int): Maximum number of samples.
        inflate (float): Inflation factor $\geq 1$ to multiply by the variance estimate to make it more conservative.
        alpha (np.ndarray): Uncertainty level in $(0,1)$. 
        control_variates (list): Integrands to use as control variates, each with the same underlying discrete distribution instance.
        control_variate_means (np.ndarray): Means of each control variate. 
    """
    self.parameters = ['abs_tol','rel_tol','n_init','n_limit','inflate','alpha','kurtmax']
    # Set Attributes
    self.abs_tol = float(abs_tol)
    self.rel_tol = float(rel_tol)
    self.n_init = float(n_init)
    self.n_limit = float(n_limit)
    self.alpha = float(alpha)
    self.inflate = float(inflate)
    self.alpha_sigma = float(self.alpha) / 2.  # the uncertainty for variance estimation
    self.kurtmax = (n_init - 3) / (n_init - 1) + \
        (self.alpha_sigma * n_init) / (1 - self.alpha_sigma) * \
        (1 - 1. / self.inflate**2)**2
    # QMCPy Objs
    self.integrand = integrand
    self.true_measure = self.integrand.true_measure
    self.discrete_distrib = self.integrand.discrete_distrib
    super(CubMCG,self).__init__(allowed_distribs=[AbstractIIDDiscreteDistribution],allow_vectorized_integrals=False)
    assert self.integrand.d_indv==()
    # control variates
    self.cv_mu = np.atleast_1d(control_variate_means)
    self.cv = control_variates
    if isinstance(self.cv,AbstractIntegrand):
        self.cv = [self.cv]
        self.cv_mu = self.cv_mu[None,...]
    assert isinstance(self.cv,list), "cv must be a list of AbstractIntegrand objects"
    for cv in self.cv:
        if (not isinstance(cv,AbstractIntegrand)) or (cv.discrete_distrib!=self.discrete_distrib) or (cv.d_indv!=self.integrand.d_indv):
            raise ParameterError('''
                    Each control variates discrete distribution must be an AbstractIntegrand instance 
                    with the same discrete distribution as the main integrand. d_indv must also match 
                    that of the main integrand instance for each control variate.''')
    self.ncv = len(self.cv)
    if self.ncv>0:
        assert self.cv_mu.shape==((self.ncv,)+self.integrand.d_indv), "Control variate means should have shape (len(control variates),d_indv)."
        self.parameters += ['cv','cv_mu']

CubMCCLT

Bases: AbstractStoppingCriterion

IID Monte Carlo stopping criterion based on the Central Limit Theorem in a two step method.

Examples:

>>> ao = FinancialOption(IIDStdUniform(52,seed=7))
>>> sc = CubMCCLT(ao,abs_tol=.05)
>>> solution,data = sc.integrate()
>>> data
Data (Data)
    solution        1.777
    bound_low       1.723
    bound_high      1.831
    bound_diff      0.107
    n_total         74235
    time_integrate  ...
CubMCCLT (AbstractStoppingCriterion)
    abs_tol         0.050
    rel_tol         0
    n_init          2^(10)
    n_limit         2^(30)
    inflate         1.200
    alpha           0.010
FinancialOption (AbstractIntegrand)
    option          ASIAN
    call_put        CALL
    volatility      2^(-1)
    start_price     30
    strike_price    35
    interest_rate   0
    t_final         1
    asian_mean      ARITHMETIC
BrownianMotion (AbstractTrueMeasure)
    time_vec        [0.019 0.038 0.058 ... 0.962 0.981 1.   ]
    drift           0
    mean            [0. 0. 0. ... 0. 0. 0.]
    covariance      [[0.019 0.019 0.019 ... 0.019 0.019 0.019]
                     [0.019 0.038 0.038 ... 0.038 0.038 0.038]
                     [0.019 0.038 0.058 ... 0.058 0.058 0.058]
                     ...
                     [0.019 0.038 0.058 ... 0.962 0.962 0.962]
                     [0.019 0.038 0.058 ... 0.962 0.981 0.981]
                     [0.019 0.038 0.058 ... 0.962 0.981 1.   ]]
    decomp_type     PCA
IIDStdUniform (AbstractIIDDiscreteDistribution)
    d               52
    replications    1
    entropy         7

Control variates 

>>> iid = IIDStdUniform(52,seed=7)
>>> ao = FinancialOption(iid,option="ASIAN")
>>> eo = FinancialOption(iid,option="EUROPEAN")
>>> lin0 = Linear0(iid)
>>> sc = CubMCCLT(ao,abs_tol=.05,
...     control_variates = [eo,lin0],
...     control_variate_means = [eo.get_exact_value(),0])
>>> solution,data = sc.integrate()
>>> data
Data (Data)
    solution        1.790
    bound_low       1.738
    bound_high      1.843
    bound_diff      0.104
    n_total         27777
    time_integrate  ...
CubMCCLT (AbstractStoppingCriterion)
    abs_tol         0.050
    rel_tol         0
    n_init          2^(10)
    n_limit         2^(30)
    inflate         1.200
    alpha           0.010
    cv              [FinancialOption (AbstractIntegrand)
                         option          EUROPEAN
                         call_put        CALL
                         volatility      2^(-1)
                         start_price     30
                         strike_price    35
                         interest_rate   0
                         t_final         1               Linear0 (AbstractIntegrand)]
    cv_mu           [4.211 0.   ]
FinancialOption (AbstractIntegrand)
    option          ASIAN
    call_put        CALL
    volatility      2^(-1)
    start_price     30
    strike_price    35
    interest_rate   0
    t_final         1
    asian_mean      ARITHMETIC
BrownianMotion (AbstractTrueMeasure)
    time_vec        [0.019 0.038 0.058 ... 0.962 0.981 1.   ]
    drift           0
    mean            [0. 0. 0. ... 0. 0. 0.]
    covariance      [[0.019 0.019 0.019 ... 0.019 0.019 0.019]
                     [0.019 0.038 0.038 ... 0.038 0.038 0.038]
                     [0.019 0.038 0.058 ... 0.058 0.058 0.058]
                     ...
                     [0.019 0.038 0.058 ... 0.962 0.962 0.962]
                     [0.019 0.038 0.058 ... 0.962 0.981 0.981]
                     [0.019 0.038 0.058 ... 0.962 0.981 1.   ]]
    decomp_type     PCA
IIDStdUniform (AbstractIIDDiscreteDistribution)
    d               52
    replications    1
    entropy         7

Parameters:

Name Type Description Default
integrand AbstractIntegrand

The integrand.

 required
abs_tol ndarray

Absolute error tolerance.

 0.01
rel_tol ndarray

Relative error tolerance.

 0.0
n_init int

Initial number of samples.

1024
n_limit int

Maximum number of samples.

 2 ** 30
inflate float

Inflation factor \(\geq 1\) to multiply by the variance estimate to make it more conservative.

 1.2
alpha ndarray

Uncertainty level in \((0,1)\).

0.01
control_variates list

Integrands to use as control variates, each with the same underlying discrete distribution instance.

 []
control_variate_means ndarray

Means of each control variate.

 []
Source code in qmcpy/stopping_criterion/cub_mc_clt.py 
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
def __init__(self, 
             integrand, 
             abs_tol = 1e-2, 
             rel_tol = 0., 
             n_init = 1024, 
             n_limit = 2**30,
             inflate = 1.2, 
             alpha = 0.01, 
             control_variates = [], 
             control_variate_means = [],
             ):
    r"""
    Args:
        integrand (AbstractIntegrand): The integrand.
        abs_tol (np.ndarray): Absolute error tolerance.
        rel_tol (np.ndarray): Relative error tolerance.
        n_init (int): Initial number of samples. 
        n_limit (int): Maximum number of samples.
        inflate (float): Inflation factor $\geq 1$ to multiply by the variance estimate to make it more conservative.
        alpha (np.ndarray): Uncertainty level in $(0,1)$. 
        control_variates (list): Integrands to use as control variates, each with the same underlying discrete distribution instance.
        control_variate_means (np.ndarray): Means of each control variate. 
    """
    self.parameters = ['abs_tol','rel_tol','n_init','n_limit','inflate','alpha']
    # Set Attributes
    self.abs_tol = abs_tol
    self.rel_tol = rel_tol
    self.n_init = n_init
    self.n_limit = n_limit
    assert self.n_limit>(2*self.n_init), "require n_limit is at least twic as much as n_init"
    self.alpha = alpha
    self.inflate = inflate
    assert self.inflate>=1
    assert 0<self.alpha<1
    # QMCPy Objs
    self.integrand = integrand
    self.true_measure = self.integrand.true_measure
    self.discrete_distrib = self.integrand.discrete_distrib
    super(CubMCCLT,self).__init__(allowed_distribs=[AbstractIIDDiscreteDistribution],allow_vectorized_integrals=True)
    assert self.integrand.d_indv==()
    # control variates
    self.cv_mu = np.atleast_1d(control_variate_means)
    self.cv = control_variates
    if isinstance(self.cv,AbstractIntegrand):
        self.cv = [self.cv]
        self.cv_mu = self.cv_mu[None,...]
    assert isinstance(self.cv,list), "cv must be a list of AbstractIntegrand objects"
    for cv in self.cv:
        if (not isinstance(cv,AbstractIntegrand)) or (cv.discrete_distrib!=self.discrete_distrib) or (cv.d_indv!=self.integrand.d_indv):
            raise ParameterError('''
                    Each control variates discrete distribution must be an AbstractIntegrand instance 
                    with the same discrete distribution as the main integrand. d_indv must also match 
                    that of the main integrand instance for each control variate.''')
    self.ncv = len(self.cv)
    if self.ncv>0:
        assert self.cv_mu.shape==((self.ncv,)+self.integrand.d_indv), "Control variate means should have shape (len(control variates),d_indv)."
        self.parameters += ['cv','cv_mu']
    self.z_star = -norm.ppf(self.alpha/2.)

Multilevel QMC

CubMLQMCCont

Bases: AbstractCubMLQMC

Multilevel Quasi-Monte Carlo stopping criterion with continuation.

Examples:

>>> fo = FinancialOption(DigitalNetB2(seed=7,replications=32))
>>> sc = CubMLQMCCont(fo,abs_tol=1e-3)
>>> solution,data = sc.integrate()
>>> data
Data (Data)
    solution        1.784
    n_total         4718592
    levels          2^(2)
    n_level         [65536 32768 32768 16384]
    mean_level      [1.718 0.051 0.012 0.003]
    var_level       [1.169e-08 2.569e-08 1.850e-08 5.209e-08]
    bias_estimate   2.78e-04
    time_integrate  ...
CubMLQMCCont (AbstractStoppingCriterion)
    rmse_tol        3.88e-04
    n_init          2^(8)
    n_limit         10000000000
    replications    2^(5)
    levels_min      2^(1)
    levels_max      10
    n_tols          10
    inflate         1.668
    theta_init      2^(-1)
    theta           2^(-3)
FinancialOption (AbstractIntegrand)
    option          ASIAN
    call_put        CALL
    volatility      2^(-1)
    start_price     30
    strike_price    35
    interest_rate   0
    t_final         1
    asian_mean      ARITHMETIC
BrownianMotion (AbstractTrueMeasure)
    time_vec        1
    drift           0
    mean            0
    covariance      1
    decomp_type     PCA
DigitalNetB2 (AbstractLDDiscreteDistribution)
    d               1
    replications    2^(5)
    randomize       LMS DS
    gen_mats_source joe_kuo.6.21201.txt
    order           RADICAL INVERSE
    t               63
    alpha           1
    n_limit         2^(32)
    entropy         7

References:

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

Parameters:

Name Type Description Default
integrand AbstractIntegrand

The integrand.

 required
abs_tol ndarray

Absolute error tolerance.

 0.05
rmse_tol ndarray

Root mean squared error tolerance. If supplied, then absolute tolerance and alpha are ignored in favor of the rmse tolerance.

None
n_init int

Initial number of samples.

256
n_limit int

Maximum number of samples.

 10000000000.0
inflate float

Coarser tolerance multiplication factor \(\geq 1\).

 100 ** (1 / 9)
alpha ndarray

Uncertainty level in \((0,1)\).

0.01
levels_min int

Minimum level of refinement \(\geq 2\).

 2
levels_max int

Maximum level of refinement \(\geq\) levels_min.

 10
n_tols int

Number of coarser tolerances to run.

 10
theta_init float

Initial error splitting constant.

 0.5
Source code in qmcpy/stopping_criterion/cub_mlqmc_cont.py 
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
def __init__(self, 
             integrand, 
             abs_tol = .05, 
             rmse_tol = None,
             n_init = 256,
             n_limit = 1e10,
             inflate = 100**(1/9),
             alpha = .01,   
             levels_min = 2, 
             levels_max = 10, 
             n_tols = 10, 
             theta_init = 0.5,
             ):
    r"""
    Args:
        integrand (AbstractIntegrand): The integrand.
        abs_tol (np.ndarray): Absolute error tolerance.
        rmse_tol (np.ndarray): Root mean squared error tolerance. 
            If supplied, then absolute tolerance and alpha are ignored in favor of the rmse tolerance. 
        n_init (int): Initial number of samples. 
        n_limit (int): Maximum number of samples.
        inflate (float): Coarser tolerance multiplication factor $\geq 1$.
        alpha (np.ndarray): Uncertainty level in $(0,1)$. 
        levels_min (int): Minimum level of refinement $\geq 2$.
        levels_max (int): Maximum level of refinement $\geq$ `levels_min`.
        n_tols (int): Number of coarser tolerances to run.
        theta_init (float): Initial error splitting constant.
    """
    self.parameters = ['rmse_tol','n_init','n_limit','replications','levels_min',
        'levels_max','n_tols','inflate','theta_init','theta']
    # initialization
    if rmse_tol:
        self.target_tol = float(rmse_tol)
    else: # use absolute tolerance
        self.target_tol =  float(abs_tol) / norm.ppf(1-alpha/2)
    self.n_init = n_init
    self.n_limit = n_limit
    self.levels_min = levels_min
    self.levels_max = levels_max
    self.theta_init = theta_init
    self.theta = theta_init
    self.n_tols = n_tols
    self.alpha = alpha
    self.inflate = inflate
    assert self.inflate>=1
    assert 0<self.alpha<1
    # QMCPy Objs
    self.integrand = integrand
    self.true_measure = self.integrand.true_measure
    self.discrete_distrib = self.integrand.discrete_distrib
    super(CubMLQMCCont,self).__init__(allowed_distribs=[AbstractLDDiscreteDistribution],allow_vectorized_integrals=False)
    self.replications = self.discrete_distrib.replications 
    assert self.replications>=4, "require at least 4 replications"

CubMLQMC

Bases: AbstractCubMLQMC

Multilevel Quasi-Monte Carlo stopping criterion.

Examples:

>>> fo = FinancialOption(DigitalNetB2(seed=7,replications=32))
>>> sc = CubMLQMC(fo,abs_tol=3e-3)
>>> solution,data = sc.integrate()
>>> data
Data (Data)
    solution        1.784
    n_total         2359296
    levels          2^(2)
    n_level         [32768 16384 16384  8192]
    mean_level      [1.718 0.051 0.012 0.003]
    var_level       [7.119e-08 1.409e-07 9.668e-08 1.852e-07]
    bias_estimate   2.99e-04
    time_integrate  ...
CubMLQMC (AbstractStoppingCriterion)
    rmse_tol        0.001
    n_init          2^(8)
    n_limit         10000000000
    replications    2^(5)
FinancialOption (AbstractIntegrand)
    option          ASIAN
    call_put        CALL
    volatility      2^(-1)
    start_price     30
    strike_price    35
    interest_rate   0
    t_final         1
    asian_mean      ARITHMETIC
BrownianMotion (AbstractTrueMeasure)
    time_vec        1
    drift           0
    mean            0
    covariance      1
    decomp_type     PCA
DigitalNetB2 (AbstractLDDiscreteDistribution)
    d               1
    replications    2^(5)
    randomize       LMS DS
    gen_mats_source joe_kuo.6.21201.txt
    order           RADICAL INVERSE
    t               63
    alpha           1
    n_limit         2^(32)
    entropy         7

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.

Parameters:

Name Type Description Default
integrand AbstractIntegrand

The integrand.

 required
abs_tol ndarray

Absolute error tolerance.

 0.05
rmse_tol ndarray

Root mean squared error tolerance. If supplied, then absolute tolerance and alpha are ignored in favor of the rmse tolerance.

None
n_init int

Initial number of samples.

256
n_limit int

Maximum number of samples.

 10000000000.0
alpha ndarray

Uncertainty level in \((0,1)\).

0.01
levels_min int

Minimum level of refinement \(\geq 2\).

 2
levels_max int

Maximum level of refinement \(\geq\) levels_min.

 10
Source code in qmcpy/stopping_criterion/cub_mlqmc.py 
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
def __init__(self, 
             integrand, 
             abs_tol = .05, 
             rmse_tol = None, 
             n_init = 256, 
             n_limit = 1e10,
             alpha = .01, 
             levels_min = 2,
             levels_max = 10,
             ):
    r"""
    Args:
        integrand (AbstractIntegrand): The integrand.
        abs_tol (np.ndarray): Absolute error tolerance.
        rmse_tol (np.ndarray): Root mean squared error tolerance. 
            If supplied, then absolute tolerance and alpha are ignored in favor of the rmse tolerance. 
        n_init (int): Initial number of samples. 
        n_limit (int): Maximum number of samples.
        alpha (np.ndarray): Uncertainty level in $(0,1)$. 
        levels_min (int): Minimum level of refinement $\geq 2$.
        levels_max (int): Maximum level of refinement $\geq$ `levels_min`.
    """
    self.parameters = ['rmse_tol','n_init','n_limit','replications']
    # initialization
    if rmse_tol:
        self.rmse_tol = float(rmse_tol)
    else: # use absolute tolerance
        self.rmse_tol =  float(abs_tol) / norm.ppf(1-alpha/2)
    self.alpha = alpha 
    assert 0<self.alpha<1
    self.n_init = n_init
    self.n_limit = n_limit
    self.levels_min = levels_min
    self.levels_max = levels_max
    # QMCPy Objs
    self.integrand = integrand
    self.true_measure = self.integrand.true_measure
    self.discrete_distrib = self.integrand.discrete_distrib
    super(CubMLQMC,self).__init__(allowed_distribs=[AbstractLDDiscreteDistribution],allow_vectorized_integrals=False)
    self.replications = self.discrete_distrib.replications 
    assert self.replications>=4, "require at least 4 replications"

Multilevel IID MC

CubMLMCCont

Bases: AbstractCubMLMC

Multilevel IID Monte Carlo stopping criterion with continuation.

Examples:

>>> fo = FinancialOption(IIDStdUniform(seed=7))
>>> sc = CubMLMCCont(fo,abs_tol=1.5e-2)
>>> solution,data = sc.integrate()
>>> data
Data (Data)
    solution        1.771
    n_total         2291120
    levels          3
    n_level         [1094715  222428   79666     912     256]
    mean_level      [1.71  0.048 0.012]
    var_level       [21.826  1.768  0.453]
    cost_per_sample [2. 4. 8.]
    alpha           1.970
    beta            1.965
    gamma           1.000
    time_integrate  ...
CubMLMCCont (AbstractStoppingCriterion)
    rmse_tol        0.006
    n_init          2^(8)
    levels_min      2^(1)
    levels_max      10
    n_tols          10
    inflate         1.668
    theta_init      2^(-1)
    theta           0.010
FinancialOption (AbstractIntegrand)
    option          ASIAN
    call_put        CALL
    volatility      2^(-1)
    start_price     30
    strike_price    35
    interest_rate   0
    t_final         1
    asian_mean      ARITHMETIC
BrownianMotion (AbstractTrueMeasure)
    time_vec        1
    drift           0
    mean            0
    covariance      1
    decomp_type     PCA
IIDStdUniform (AbstractIIDDiscreteDistribution)
    d               1
    replications    1
    entropy         7

References:

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

Parameters:

Name Type Description Default
integrand AbstractIntegrand

The integrand.

 required
abs_tol ndarray

Absolute error tolerance.

 0.05
rmse_tol ndarray

Root mean squared error tolerance. If supplied, then absolute tolerance and alpha are ignored in favor of the rmse tolerance.

None
n_init int

Initial number of samples.

256
n_limit int

Maximum number of samples.

 10000000000.0
inflate float

Coarser tolerance multiplication factor \(\geq 1\).

 100 ** (1 / 9)
alpha ndarray

Uncertainty level in \((0,1)\).

0.01
levels_min int

Minimum level of refinement \(\geq 2\).

 2
levels_max int

Maximum level of refinement \(\geq\) levels_min.

 10
n_tols int

Number of coarser tolerances to run.

 10
theta_init float

Initial error splitting constant.

 0.5
Source code in qmcpy/stopping_criterion/cub_mlmc_cont.py 
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
def __init__(self, 
             integrand, 
             abs_tol = .05, 
             rmse_tol = None, 
             n_init = 256, 
             n_limit = 1e10,
             inflate = 100**(1/9),
             alpha = .01, 
             levels_min = 2,
             levels_max = 10, 
             n_tols = 10, 
             theta_init = 0.5,
             ):
    r"""
    Args:
        integrand (AbstractIntegrand): The integrand.
        abs_tol (np.ndarray): Absolute error tolerance.
        rmse_tol (np.ndarray): Root mean squared error tolerance. 
            If supplied, then absolute tolerance and alpha are ignored in favor of the rmse tolerance. 
        n_init (int): Initial number of samples. 
        n_limit (int): Maximum number of samples.
        inflate (float): Coarser tolerance multiplication factor $\geq 1$.
        alpha (np.ndarray): Uncertainty level in $(0,1)$. 
        levels_min (int): Minimum level of refinement $\geq 2$.
        levels_max (int): Maximum level of refinement $\geq$ `levels_min`.
        n_tols (int): Number of coarser tolerances to run.
        theta_init (float): Initial error splitting constant.
    """
    self.parameters = ['rmse_tol','n_init','levels_min','levels_max','n_tols',
        'inflate','theta_init','theta']
    if levels_min < 2:
        raise ParameterError('needs levels_min >= 2')
    if levels_max < levels_min:
        raise ParameterError('needs levels_max >= levels_min')
    if n_init <= 0:
        raise ParameterError('needs n_init > 0')
    # initialization
    if rmse_tol:
        self.target_tol = float(rmse_tol)
    else: # use absolute tolerance
        self.target_tol =  float(abs_tol) / norm.ppf(1-alpha/2)
    self.n_init = n_init
    self.n_limit = n_limit
    self.levels_min = levels_min
    self.levels_max = levels_max
    self.theta_init = theta_init
    self.theta = theta_init
    self.n_tols = n_tols
    self.inflate = inflate
    self.integrand = integrand
    self.true_measure = self.integrand.true_measure
    self.discrete_distrib = self.integrand.discrete_distrib
    self.alpha0 = -1 
    self.beta0 = -1 
    self.gamma0 = -1
    self.alpha = alpha
    self.inflate = inflate
    assert self.inflate>=1
    assert 0<self.alpha<1
    super(CubMLMCCont,self).__init__(allowed_distribs=[AbstractIIDDiscreteDistribution],allow_vectorized_integrals=False)

CubMLMC

Bases: AbstractCubMLMC

Multilevel IID Monte Carlo stopping criterion.

Examples:

>>> fo = FinancialOption(IIDStdUniform(seed=7))
>>> sc = CubMLMC(fo,abs_tol=1.5e-2)
>>> solution,data = sc.integrate()
>>> data
Data (Data)
    solution        1.785
    n_total         3577556
    levels          2^(2)
    n_level         [2438191  490331  207606   62905]
    mean_level      [1.715 0.053 0.013 0.003]
    var_level       [21.829  1.761  0.452  0.11 ]
    cost_per_sample [ 2.  4.  8. 16.]
    alpha           2.008
    beta            1.997
    gamma           1.000
    time_integrate  ...
CubMLMC (AbstractStoppingCriterion)
    rmse_tol        0.006
    n_init          2^(8)
    levels_min      2^(1)
    levels_max      10
    theta           2^(-1)
FinancialOption (AbstractIntegrand)
    option          ASIAN
    call_put        CALL
    volatility      2^(-1)
    start_price     30
    strike_price    35
    interest_rate   0
    t_final         1
    asian_mean      ARITHMETIC
BrownianMotion (AbstractTrueMeasure)
    time_vec        1
    drift           0
    mean            0
    covariance      1
    decomp_type     PCA
IIDStdUniform (AbstractIIDDiscreteDistribution)
    d               1
    replications    1
    entropy         7

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.

  2. http://people.maths.ox.ac.uk/~gilesm/mlmc/#MATLAB.

Parameters:

Name Type Description Default
integrand AbstractIntegrand

The integrand.

 required
abs_tol ndarray

Absolute error tolerance.

 0.05
rmse_tol ndarray

Root mean squared error tolerance. If supplied, then absolute tolerance and alpha are ignored in favor of the rmse tolerance.

None
n_init int

Initial number of samples.

256
n_limit int

Maximum number of samples.

 10000000000.0
alpha ndarray

Uncertainty level in \((0,1)\).

0.01
levels_min int

Minimum level of refinement \(\geq 2\).

 2
levels_max int

Maximum level of refinement \(\geq\) levels_min.

 10
alpha0 float

Weak error is \(\mathcal{O}(2^{-\alpha_0\ell})\) in the level \(\ell\). If alpha0\(\leq 0\) then it will be estimated.

-1.0
beta0 float

Variance is \(\mathcal{O}(2^{-\beta_0\ell})\) in the level \(\ell\). If beta0\(\leq 0\) then it will be estimated.

-1.0
gamma0 float

Sample cost is \(\mathcal{O}(2^{\gamma_0\ell})\) in the level \(\ell\). If gamma0\(\leq 0\) then it will be estimated.

 -1.0
Source code in qmcpy/stopping_criterion/cub_mlmc.py 
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
def __init__(self, 
             integrand, 
             abs_tol = .05, 
             rmse_tol = None, 
             n_init = 256, 
             n_limit = 1e10, 
             alpha = .01, 
             levels_min = 2, 
             levels_max = 10, 
             alpha0 = -1., 
             beta0 = -1., 
             gamma0 = -1.,
             ):
    r"""
    Args:
        integrand (AbstractIntegrand): The integrand.
        abs_tol (np.ndarray): Absolute error tolerance.
        rmse_tol (np.ndarray): Root mean squared error tolerance. 
            If supplied, then absolute tolerance and alpha are ignored in favor of the rmse tolerance. 
        n_init (int): Initial number of samples. 
        n_limit (int): Maximum number of samples.
        alpha (np.ndarray): Uncertainty level in $(0,1)$. 
        levels_min (int): Minimum level of refinement $\geq 2$.
        levels_max (int): Maximum level of refinement $\geq$ `levels_min`.
        alpha0 (float): Weak error is $\mathcal{O}(2^{-\alpha_0\ell})$ in the level $\ell$. If `alpha0`$\leq 0$ then it will be estimated. 
        beta0 (float): Variance is $\mathcal{O}(2^{-\beta_0\ell})$ in the level $\ell$. If `beta0`$\leq 0$ then it will be estimated. 
        gamma0 (float): Sample cost is $\mathcal{O}(2^{\gamma_0\ell})$ in the level $\ell$. If `gamma0`$\leq 0$ then it will be estimated. 
    """
    self.parameters = ['rmse_tol','n_init','levels_min','levels_max','theta']
    if levels_min < 2:
        raise ParameterError('needs levels_min >= 2')
    if levels_max < levels_min:
        raise ParameterError('needs levels_max >= levels_min')
    if n_init <= 0:
        raise ParameterError('needs n_init > 0')
    # initialization
    if rmse_tol:
        self.rmse_tol = float(rmse_tol)
    else: # use absolute tolerance
        self.rmse_tol =  float(abs_tol) / norm.ppf(1-alpha/2)
    self.alpha = alpha 
    assert 0<self.alpha<1
    self.n_init = n_init
    self.n_limit = n_limit
    self.levels_min = levels_min
    self.levels_max = levels_max
    self.theta = 0.5
    self.alpha0 = alpha0
    self.beta0 = beta0
    self.gamma0 = gamma0
    # QMCPy Objs
    self.integrand = integrand
    self.true_measure = self.integrand.true_measure
    self.discrete_distrib = self.integrand.discrete_distrib
    super(CubMLMC,self).__init__(allowed_distribs=[AbstractIIDDiscreteDistribution], allow_vectorized_integrals=False)

Failure Probability

PFGPCI

Bases: AbstractStoppingCriterion

Probability of failure estimation using adaptive Gaussian process construction and resulting credible intervals.

Examples:

>>> pfgpci = PFGPCI(
...     integrand = Ishigami(DigitalNetB2(3,seed=17)),
...     failure_threshold = 0,
...     failure_above_threshold = False, 
...     abs_tol = 2.5e-2,
...     alpha = 1e-1,
...     n_init = 64,
...     init_samples = None,
...     batch_sampler = PFSampleErrorDensityAR(verbose=False),
...     n_batch = 16,
...     n_limit = 128,
...     n_approx = 2**18,
...     gpytorch_prior_mean = gpytorch.means.ZeroMean(),
...     gpytorch_prior_cov = gpytorch.kernels.ScaleKernel(
...         gpytorch.kernels.MaternKernel(nu=2.5,lengthscale_constraint = gpytorch.constraints.Interval(.5,1)),
...         outputscale_constraint = gpytorch.constraints.Interval(1e-8,.5)),
...     gpytorch_likelihood = gpytorch.likelihoods.GaussianLikelihood(noise_constraint = gpytorch.constraints.Interval(1e-12,1e-8)),
...     gpytorch_marginal_log_likelihood_func = lambda likelihood,gpyt_model: gpytorch.mlls.ExactMarginalLogLikelihood(likelihood,gpyt_model),
...     torch_optimizer_func = lambda gpyt_model: torch.optim.Adam(gpyt_model.parameters(),lr=0.1),
...     gpytorch_train_iter = 100,
...     gpytorch_use_gpu = False,
...     verbose = False,
...     n_ref_approx = 2**22,
...     seed_ref_approx = 11)
>>> solution,data = pfgpci.integrate(seed=7,refit=True)
>>> data
PFGPCIData (Data)
    solution        0.158
    error_bound     0.022
    bound_low       0.136
    bound_high      0.180
    n_total         112
    time_integrate  ...
PFGPCI (AbstractStoppingCriterion)
    abs_tol         0.025
    n_init          2^(6)
    n_limit         2^(7)
    n_batch         2^(4)
Ishigami (AbstractIntegrand)
Uniform (AbstractTrueMeasure)
    lower_bound     -3.142
    upper_bound     3.142
DigitalNetB2 (AbstractLDDiscreteDistribution)
    d               3
    replications    1
    randomize       LMS DS
    gen_mats_source joe_kuo.6.21201.txt
    order           RADICAL INVERSE
    t               63
    alpha           1
    n_limit         2^(32)
    entropy         17
>>> df = data.get_results_dict()
>>> with np.printoptions(formatter={"float": lambda x: "%-10.2e"%x, "int": lambda x: "%-10d"%x, "bool": lambda x: "%-10s"%x}):
...     for k,v in df.items():
...         print("%15s: %s"%(k,str(v)))
           iter: [0          1          2          3         ]
          n_sum: [64         80         96         112       ]
        n_batch: [64         16         16         16        ]
   error_bounds: [5.58e-02   3.92e-02   3.05e-02   2.16e-02  ]
         ci_low: [8.66e-02   1.16e-01   1.19e-01   1.36e-01  ]
        ci_high: [1.98e-01   1.95e-01   1.80e-01   1.80e-01  ]
      solutions: [1.42e-01   1.55e-01   1.50e-01   1.58e-01  ]
  solutions_ref: [1.62e-01   1.62e-01   1.62e-01   1.62e-01  ]
      error_ref: [2.01e-02   7.02e-03   1.28e-02   4.52e-03  ]
          in_ci: [True       True       True       True      ]

References:

  1. Sorokin, Aleksei G., and Vishwas Rao.
    "Credible Intervals for Probability of Failure with Gaussian Processes."
    arXiv preprint arXiv:2311.07733 (2023).

Parameters:

Name Type Description Default
integrand AbstractIntegrand

The integrand.

 required
failure_threshold float

Thresholds for failure.

required
failure_above_threshold bool

Set to True if failure occurs when the simulation exceeds failure_threshold and False otherwise.

 required
abs_tol float

The desired maximum distance from the estimate to either end of the credible interval.

 0.005
n_init float

Initial number of samples from integrand.discrete_distrib from which to build the first surrogate GP

 64
n_limit int

Budget of simulations.

 1000
n_batch int

The number of samples per batch to draw from batch_sampler.

4
alpha float

The credible interval is constructed to hold with probability at least 1 - alpha

 0.01
init_samples float

If the simulation has already been run, pass in (x,y) where x are past samples from the discrete distribution and y are corresponding simulation evaluations.

None
batch_sampler Suggester or AbstractDiscreteDistribution

A suggestion scheme for future samples.

PFSampleErrorDensityAR()
n_approx int

Number of points from integrand.discrete_distrib used to approximate estimate and credible interval bounds

 2 ** 20
gpytorch_prior_mean means

prior mean function of the GP

 ZeroMean()
gpytorch_prior_cov kernels

Prior covariance kernel of the GP

 ScaleKernel(MaternKernel(nu=2.5))
gpytorch_likelihood likelihoods

GP likelihood, require one of gpytorch.likelihoods.{GaussianLikelihood, GaussianLikelihoodWithMissingObs, FixedNoiseGaussianLikelihood}

 GaussianLikelihood(noise_constraint=Interval(1e-12, 1e-08))
gpytorch_marginal_log_likelihood_func callable

Function taking in the likelihood and gpytorch model and returning a marginal log likelihood from gpytorch.mlls

 lambda likelihood, gpyt_model: ExactMarginalLogLikelihood(likelihood, gpyt_model)
torch_optimizer_func callable

Function taking in the gpytorch model and returning an optimizer from torch.optim

 lambda gpyt_model: Adam(parameters(), lr=0.1)
gpytorch_train_iter int

Training iterations for the GP in gpytorch

 100
gpytorch_use_gpu bool

If True, have gpytorch use a GPU for fitting and trining the GP

 False
verbose int

If verbose > 0, print information through the call to integrate()

 False
n_ref_approx int

If n_ref_approx > 0, use n_ref_approx points to get a reference QMC approximation of the true solution. Caution: If n_ref_approx > 0, it should be a large int e.g. 2**22, in which case it is only helpful for cheap to evaluate simulations

2 ** 22
seed_ref_approx int

Seed for the reference approximation. Only applies when n_ref_approx>0

 None
Source code in qmcpy/stopping_criterion/pf_gp_ci.py 
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
def __init__(self, 
    integrand,
    failure_threshold,
    failure_above_threshold,
    abs_tol = 5e-3,
    n_init = 64,
    n_limit = 1000,
    alpha = 1e-2,
    init_samples = None,
    batch_sampler = PFSampleErrorDensityAR(),
    n_batch = 4,
    n_approx = 2**20,
    gpytorch_prior_mean = gpytorch.means.ZeroMean(),
    gpytorch_prior_cov = gpytorch.kernels.ScaleKernel(gpytorch.kernels.MaternKernel(nu = 2.5)),
    gpytorch_likelihood = gpytorch.likelihoods.GaussianLikelihood(noise_constraint = gpytorch.constraints.Interval(1e-12,1e-8)),
    gpytorch_marginal_log_likelihood_func = lambda likelihood,gpyt_model: gpytorch.mlls.ExactMarginalLogLikelihood(likelihood,gpyt_model),
    torch_optimizer_func = lambda gpyt_model: torch.optim.Adam(gpyt_model.parameters(),lr=0.1),
    gpytorch_train_iter = 100,
    gpytorch_use_gpu = False,
    verbose = False,
    n_ref_approx = 2**22,
    seed_ref_approx = None):
    '''
    Args:
        integrand (AbstractIntegrand): The integrand.
        failure_threshold (float): Thresholds for failure. 
        failure_above_threshold (bool): Set to `True` if failure occurs when the simulation exceeds `failure_threshold` and False otherwise.
        abs_tol (float): The desired maximum distance from the estimate to either end of the credible interval.
        n_init (float): Initial number of samples from integrand.discrete_distrib from which to build the first surrogate GP
        n_limit (int): Budget of simulations.
        n_batch (int): The number of samples per batch to draw from batch_sampler. 
        alpha (float): The credible interval is constructed to hold with probability at least 1 - alpha
        init_samples (float): If the simulation has already been run, pass in (x,y) where x are past samples from the discrete distribution and y are corresponding simulation evaluations. 
        batch_sampler (Suggester or AbstractDiscreteDistribution): A suggestion scheme for future samples. 
        n_approx (int): Number of points from integrand.discrete_distrib used to approximate estimate and credible interval bounds
        gpytorch_prior_mean (gpytorch.means): prior mean function of the GP
        gpytorch_prior_cov (gpytorch.kernels): Prior covariance kernel of the GP
        gpytorch_likelihood (gpytorch.likelihoods): GP likelihood, require one of gpytorch.likelihoods.{GaussianLikelihood, GaussianLikelihoodWithMissingObs, FixedNoiseGaussianLikelihood}
        gpytorch_marginal_log_likelihood_func (callable): Function taking in the likelihood and gpytorch model and returning a marginal log likelihood from gpytorch.mlls
        torch_optimizer_func (callable): Function taking in the gpytorch model and returning an optimizer from torch.optim
        gpytorch_train_iter (int): Training iterations for the GP in gpytorch
        gpytorch_use_gpu (bool): If True, have gpytorch use a GPU for fitting and trining the GP
        verbose (int): If verbose > 0, print information through the call to integrate()
        n_ref_approx (int): If n_ref_approx > 0, use n_ref_approx points to get a reference QMC approximation of the true solution. 
            Caution: If n_ref_approx > 0, it should be a large int e.g. 2**22, in which case it is only helpful for cheap to evaluate simulations 
        seed_ref_approx (int): Seed for the reference approximation. Only applies when n_ref_approx>0
    '''
    self.parameters = ['abs_tol','n_init','n_limit','n_batch']
    self.integrand = integrand
    self.true_measure = self.integrand.true_measure
    self.discrete_distrib = self.integrand.discrete_distrib
    self.d = self.integrand.d
    self.sampler = self.d
    self.failure_threshold = failure_threshold
    self.failure_above_threshold = failure_above_threshold
    self.abs_tol = abs_tol
    self.alpha = alpha
    assert 0<self.alpha<1
    self.n_init = n_init
    self.init_samples = init_samples is not None
    if self.init_samples: 
        self.x_init,self.y_init = init_samples
        assert self.x_init.ndim==2 and self.y_init.ndim==1
        assert self.x_init.shape[1]==self.d and len(self.y_init)==len(self.x_init)
        assert self.n_init==len(self.x_init)
        self.ytf_init = self._affine_tf(self.y_init)
    self.batch_sampler = batch_sampler
    self.n_batch = n_batch
    self.n_limit = n_limit
    assert self.n_limit >= self.n_init
    self.n_approx = n_approx
    assert (self.n_approx+self.n_init) <= 2**32
    self.gpytorch_prior_mean = gpytorch_prior_mean
    self.gpytorch_prior_cov = gpytorch_prior_cov
    self.gpytorch_likelihood = gpytorch_likelihood
    self.gpytorch_marginal_log_likelihood_func = gpytorch_marginal_log_likelihood_func
    self.torch_optimizer_func = torch_optimizer_func
    self.gpytorch_train_iter = gpytorch_train_iter
    self.gpytorch_use_gpu = gpytorch_use_gpu
    self.verbose = verbose
    self.approx_true_solution = n_ref_approx>0
    if self.approx_true_solution: 
        x = DigitalNetB2(self.d,order="GRAY",seed=seed_ref_approx)(n_ref_approx)
        y = self.integrand.f(x).squeeze()
        ytf = self._affine_tf(y)
        self.ref_approx = (ytf>=0).mean(0)
        if self.verbose: print('reference approximation with d=%d: %s'%(self.d,self.ref_approx))
    super(PFGPCI,self).__init__(allowed_distribs=[AbstractDiscreteDistribution],allow_vectorized_integrals=False)

UML Specific