Stopping Criteria

UML Overview

`AbstractStoppingCriterion`

Bases: object

Initialize a stopping criterion base class.

Parameters:

Name	Type	Description	Default
`allowed_distribs`	`list`	Allowed discrete distribution classes.	required
`allow_vectorized_integrals`	`bool`	Whether to allow integrands with vectorized outputs, meaning `integrand.d_indv != ()`.	required

Source code in qmcpy/stopping_criterion/abstract_stopping_criterion.py

def __init__(self, allowed_distribs, allow_vectorized_integrals):
    """Initialize a stopping criterion base class.

    Args:
        allowed_distribs (list): Allowed discrete distribution classes.
        allow_vectorized_integrals (bool): Whether to allow integrands with
            vectorized outputs, meaning ``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 = []
    self.elapsed_time = float(getattr(self, "elapsed_time", 0.0))

integrate

integrate(resume=None)

Determine the samples needed to satisfy the target tolerance.

Parameters:

Name	Type	Description	Default
`resume`	`Data`	Existing integration state to resume from, if supported. A valid resume checkpoint must continue the same numerical experiment without duplicating samples, losing accumulated statistics, or weakening the requested tolerance guarantee. Supported resume implementations validate and copy the checkpoint before mutating any state so the caller's saved object is preserved. Defaults to None.	`None`

Returns:

Type	Description
`tuple`	tuple[Union[float, np.ndarray], Data]: Approximation to the integral with shape `integrand.d_comb` and the corresponding data object.

Source code in qmcpy/stopping_criterion/abstract_stopping_criterion.py

def integrate(self, resume=None) -> tuple:
    """Determine the samples needed to satisfy the target tolerance.

    Args:
        resume (Data, optional): Existing integration state to resume from,
            if supported. A valid resume checkpoint must continue the same
            numerical experiment without duplicating samples, losing
            accumulated statistics, or weakening the requested tolerance
            guarantee. Supported resume implementations validate and copy
            the checkpoint before mutating any state so the caller's saved
            object is preserved. Defaults to None.

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

get_iteration_log

get_iteration_log(history=None, printed_only=True, drop_empty_columns=True, formatted=True, view='all')

Return the latest iteration log as a pandas DataFrame.

Parameters:

Name	Type	Description	Default
`history`	`list[dict] \| None`	Iteration history to format. If `None`, uses `self.iteration_history` when available.	`None`
`printed_only`	`bool`	If `True`, include only rows that were selected for printed output.	`True`
`drop_empty_columns`	`bool`	If `True`, drop columns with no values.	`True`
`formatted`	`bool`	If `True`, return formatted display values when available.	`True`
`view`	`str`	Which log view to return. `"all"` and `"current"` both return the full current log, `"without_resume"` removes `RESUME` rows, and `"stage_last"` returns one stop row per stage by taking the last non-`RESUME` row available for that stage.	`'all'`

Returns:

Type	Description
`DataFrame`	pandas.DataFrame: DataFrame representation of the iteration log.

Source code in qmcpy/stopping_criterion/abstract_stopping_criterion.py

def get_iteration_log(
    self,
    history=None,
    printed_only=True,
    drop_empty_columns=True,
    formatted=True,
    view="all",
) -> "pandas.DataFrame":
    """Return the latest iteration log as a pandas DataFrame.

    Args:
        history (list[dict] | None): Iteration history to format. If ``None``,
            uses ``self.iteration_history`` when available.
        printed_only (bool): If ``True``, include only rows that were
            selected for printed output.
        drop_empty_columns (bool): If ``True``, drop columns with no values.
        formatted (bool): If ``True``, return formatted display values when
            available.
        view (str): Which log view to return. ``"all"`` and ``"current"``
            both return the full current log, ``"without_resume"`` removes
            ``RESUME`` rows, and ``"stage_last"`` returns one stop row per
            stage by taking the last non-``RESUME`` row available for that
            stage.

    Returns:
        pandas.DataFrame: DataFrame representation of the iteration log.
    """
    self._validate_iteration_log_view(view)
    use_cache = (
        history is None
        and printed_only
        and drop_empty_columns
        and formatted
        and view in ("all", "current")
    )
    if use_cache:
        cached_history_df = getattr(self, "history_df", None)
        if cached_history_df is not None:
            return cached_history_df
    if history is None:
        history = getattr(self, "iteration_history", None)
    result = _get_iteration_log_frame(
        history,
        stopping_criterion=self,
        printed_only=printed_only,
        drop_empty_columns=False,
        formatted=formatted,
    )
    result = self._apply_iteration_log_view(result, view)
    if drop_empty_columns and not result.empty:
        result = result.dropna(axis=1, how="all")
    if use_cache:
        # Cache for future calls and populate the last data object so that
        # data.history_df is available after calling get_iteration_log().
        self.history_df = result
        last_data = getattr(self, "_last_finalized_data", None)
        if last_data is not None:
            last_data.history_df = result
    return result

format_iteration_log

format_iteration_log(history=None, printed_only=True, include_header=True)

Return the iteration log as formatted text.

Parameters:

Name	Type	Description	Default
`history`	`IterationHistoryTable \| None`	Iteration history to format. If `None`, uses `self.iteration_history` when available. Defaults to None.	`None`
`printed_only`	`bool`	If `True`, include only rows that were selected for printed output. Defaults to True.	`True`
`include_header`	`bool`	If `True`, include the trace label header before the table. Defaults to True.	`True`

Returns:

Name	Type	Description
`str`	`str`	Formatted iteration log text.

Source code in qmcpy/stopping_criterion/abstract_stopping_criterion.py

def format_iteration_log(self, history=None, printed_only=True, include_header=True) -> str:
    """Return the iteration log as formatted text.

    Args:
        history (IterationHistoryTable | None, optional): Iteration history
            to format. If ``None``, uses ``self.iteration_history`` when
            available. Defaults to None.
        printed_only (bool, optional): If ``True``, include only rows that
            were selected for printed output. Defaults to True.
        include_header (bool, optional): If ``True``, include the trace
            label header before the table. Defaults to True.

    Returns:
        str: Formatted iteration log text.
    """
    if history is None:
        history = getattr(self, "iteration_history", None)
    return _format_iteration_log(
        history,
        stopping_criterion=self,
        printed_only=printed_only,
        include_header=include_header,
    )

print_iteration_log

print_iteration_log(history=None, printed_only=True, include_header=True, file=None)

Print the iteration log for the latest run or supplied history.

Parameters:

Name	Type	Description	Default
`history`	`IterationHistoryTable \| None`	Iteration history to print. If `None`, uses `self.iteration_history` when available. Defaults to None.	`None`
`printed_only`	`bool`	If `True`, print only rows that were selected for printed output. Defaults to True.	`True`
`include_header`	`bool`	If `True`, include the trace label header before the table. Defaults to True.	`True`
`file`	`TextIO \| None`	Output stream. Defaults to `sys.stdout` when None.	`None`

Returns:

Name	Type	Description
`None`	`None`	This method writes output to `file`.

Source code in qmcpy/stopping_criterion/abstract_stopping_criterion.py

def print_iteration_log(self, history=None, printed_only=True, include_header=True, file=None) -> None:
    """Print the iteration log for the latest run or supplied history.

    Args:
        history (IterationHistoryTable | None, optional): Iteration history
            to print. If ``None``, uses ``self.iteration_history`` when
            available. Defaults to None.
        printed_only (bool, optional): If ``True``, print only rows that
            were selected for printed output. Defaults to True.
        include_header (bool, optional): If ``True``, include the trace
            label header before the table. Defaults to True.
        file (typing.TextIO | None, optional): Output stream. Defaults to
            ``sys.stdout`` when None.

    Returns:
        None: This method writes output to ``file``.
    """
    if history is None:
        history = getattr(self, "iteration_history", None)
    return _print_iteration_log(
        history,
        stopping_criterion=self,
        printed_only=printed_only,
        include_header=include_header,
        file=file,
    )

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, when supported. If supplied, reset it; otherwise ignore it.	`None`
`rel_tol`	`float`	Relative tolerance, when supported. If supplied, reset it; otherwise ignore it.	`None`
`rmse_tol`	`float`	RMSE tolerance, when supported. If supplied, reset it; otherwise ignore it. 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

def set_tolerance(self, abs_tol=None, rel_tol=None, rmse_tol=None):
    """Reset the tolerances.

    Args:
        abs_tol (float): Absolute tolerance, when supported. If supplied,
            reset it; otherwise ignore it.
        rel_tol (float): Relative tolerance, when supported. If supplied,
            reset it; otherwise ignore it.
        rmse_tol (float): RMSE tolerance, when supported. If supplied,
            reset it; otherwise ignore it. If ``rmse_tol`` is not supplied
            but ``abs_tol`` is, then ``rmse_tol = abs_tol / norm.ppf(1 -
            alpha / 2)``.
    """
    raise MethodImplementationError(self, "set_tolerance")

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:

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.
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.	`_default_fudge`
`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.	`None`
`control_variate_means`	`ndarray`	Means of each control variate.	`None`
`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

def __init__(
    self,
    integrand,
    abs_tol=1e-2,
    rel_tol=0.0,
    n_init=2**10,
    n_limit=2**35,
    error_fun="EITHER",
    fudge=_default_fudge,
    check_cone=False,
    control_variates=None,
    control_variate_means=None,
    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.
    """
    if control_variates is None:
        control_variates = []
    if control_variate_means is None:
        control_variate_means = []
    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:

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.
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.	`_default_fudge`
`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

def __init__(
    self,
    integrand,
    abs_tol=1e-2,
    rel_tol=0.0,
    n_init=2**10,
    n_limit=2**30,
    error_fun="EITHER",
    fudge=_default_fudge,
    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`.
    """
    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:

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.
Jagadeeswaran Rathinavel, Fast automatic Bayesian cubature using matching kernels and designs, PhD thesis, Illinois Institute of Technology, 2019.
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

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:

Jagadeeswaran, Rathinavel, and Fred J. Hickernell. "Fast automatic Bayesian cubature using lattice sampling." Statistics and Computing 29.6 (2019): 1215-1229.
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.
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

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

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

def __init__(
    self,
    integrand,
    abs_tol=1e-2,
    rel_tol=0.0,
    n_init=256.0,
    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

def __init__(
    self,
    integrand,
    abs_tol=1e-2,
    rel_tol=0.0,
    n_init=256.0,
    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)
    self.error_fun, _ = self._resolve_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
GeometricBrownianMotion (AbstractTrueMeasure)
    time_vec        [0.019 0.038 0.058 ... 0.962 0.981 1.   ]
    drift           0
    diffusion       2^(-2)
    mean_gbm        [30. 30. 30. ... 30. 30. 30.]
    covariance_gbm  [[  4.337   4.337   4.337 ...   4.337   4.337   4.337]
                     [  4.337   8.696   8.696 ...   8.696   8.696   8.696]
                     [  4.337   8.696  13.075 ...  13.075  13.075  13.075]
                     ...
                     [  4.337   8.696  13.075 ... 244.564 244.564 244.564]
                     [  4.337   8.696  13.075 ... 244.564 250.08  250.08 ]
                     [  4.337   8.696  13.075 ... 244.564 250.08  255.623]]
    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
GeometricBrownianMotion (AbstractTrueMeasure)
    time_vec        [0.019 0.038 0.058 ... 0.962 0.981 1.   ]
    drift           0
    diffusion       2^(-2)
    mean_gbm        [30. 30. 30. ... 30. 30. 30.]
    covariance_gbm  [[  4.337   4.337   4.337 ...   4.337   4.337   4.337]
                     [  4.337   8.696   8.696 ...   8.696   8.696   8.696]
                     [  4.337   8.696  13.075 ...  13.075  13.075  13.075]
                     ...
                     [  4.337   8.696  13.075 ... 244.564 244.564 244.564]
                     [  4.337   8.696  13.075 ... 244.564 250.08  250.08 ]
                     [  4.337   8.696  13.075 ... 244.564 250.08  255.623]]
    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
GeometricBrownianMotion (AbstractTrueMeasure)
    time_vec        [0.019 0.038 0.058 ... 0.962 0.981 1.   ]
    drift           0
    diffusion       2^(-2)
    mean_gbm        [30. 30. 30. ... 30. 30. 30.]
    covariance_gbm  [[  4.337   4.337   4.337 ...   4.337   4.337   4.337]
                     [  4.337   8.696   8.696 ...   8.696   8.696   8.696]
                     [  4.337   8.696  13.075 ...  13.075  13.075  13.075]
                     ...
                     [  4.337   8.696  13.075 ... 244.564 244.564 244.564]
                     [  4.337   8.696  13.075 ... 244.564 250.08  250.08 ]
                     [  4.337   8.696  13.075 ... 244.564 250.08  255.623]]
    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
GeometricBrownianMotion (AbstractTrueMeasure)
    time_vec        [0.019 0.038 0.058 ... 0.962 0.981 1.   ]
    drift           0
    diffusion       2^(-2)
    mean_gbm        [30. 30. 30. ... 30. 30. 30.]
    covariance_gbm  [[  4.337   4.337   4.337 ...   4.337   4.337   4.337]
                     [  4.337   8.696   8.696 ...   8.696   8.696   8.696]
                     [  4.337   8.696  13.075 ...  13.075  13.075  13.075]
                     ...
                     [  4.337   8.696  13.075 ... 244.564 244.564 244.564]
                     [  4.337   8.696  13.075 ... 244.564 250.08  250.08 ]
                     [  4.337   8.696  13.075 ... 244.564 250.08  255.623]]
    decomp_type     PCA
IIDStdUniform (AbstractIIDDiscreteDistribution)
    d               52
    replications    1
    entropy         7

References:

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
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.	`None`
`control_variate_means`	`ndarray`	Means of each control variate.	`None`

Source code in qmcpy/stopping_criterion/cub_mc_g.py

def __init__(
    self,
    integrand,
    abs_tol=1e-2,
    rel_tol=0.0,
    n_init=1024,
    n_limit=2**30,
    inflate=1.2,
    alpha=0.01,
    control_variates=None,
    control_variate_means=None,
):
    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.
    """
    if control_variates is None:
        control_variates = []
    if control_variate_means is None:
        control_variate_means = []
    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.0
    )  # the uncertainty for variance estimation
    self.kurtmax = (n_init - 3) / (n_init - 1) + (self.alpha_sigma * n_init) / (
        1 - self.alpha_sigma
    ) * (1 - 1.0 / 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._init_control_variates(control_variates, control_variate_means)
    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
GeometricBrownianMotion (AbstractTrueMeasure)
    time_vec        [0.019 0.038 0.058 ... 0.962 0.981 1.   ]
    drift           0
    diffusion       2^(-2)
    mean_gbm        [30. 30. 30. ... 30. 30. 30.]
    covariance_gbm  [[  4.337   4.337   4.337 ...   4.337   4.337   4.337]
                     [  4.337   8.696   8.696 ...   8.696   8.696   8.696]
                     [  4.337   8.696  13.075 ...  13.075  13.075  13.075]
                     ...
                     [  4.337   8.696  13.075 ... 244.564 244.564 244.564]
                     [  4.337   8.696  13.075 ... 244.564 250.08  250.08 ]
                     [  4.337   8.696  13.075 ... 244.564 250.08  255.623]]
    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
GeometricBrownianMotion (AbstractTrueMeasure)
    time_vec        [0.019 0.038 0.058 ... 0.962 0.981 1.   ]
    drift           0
    diffusion       2^(-2)
    mean_gbm        [30. 30. 30. ... 30. 30. 30.]
    covariance_gbm  [[  4.337   4.337   4.337 ...   4.337   4.337   4.337]
                     [  4.337   8.696   8.696 ...   8.696   8.696   8.696]
                     [  4.337   8.696  13.075 ...  13.075  13.075  13.075]
                     ...
                     [  4.337   8.696  13.075 ... 244.564 244.564 244.564]
                     [  4.337   8.696  13.075 ... 244.564 250.08  250.08 ]
                     [  4.337   8.696  13.075 ... 244.564 250.08  255.623]]
    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.	`None`
`control_variate_means`	`ndarray`	Means of each control variate.	`None`

Source code in qmcpy/stopping_criterion/cub_mc_clt.py

def __init__(
    self,
    integrand,
    abs_tol=1e-2,
    rel_tol=0.0,
    n_init=1024,
    n_limit=2**30,
    inflate=1.2,
    alpha=0.01,
    control_variates=None,
    control_variate_means=None,
):
    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.
    """
    if control_variates is None:
        control_variates = []
    if control_variate_means is None:
        control_variate_means = []
    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._init_control_variates(control_variates, control_variate_means)
    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.0)

Multilevel QMC

`CubMLQMCCont`

Bases: AbstractCubMLQMC

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

def __init__(
    self,
    integrand,
    abs_tol=0.05,
    rmse_tol=None,
    n_init=256,
    n_limit=1e10,
    inflate=100 ** (1 / 9),
    alpha=0.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_rmse_tol = float(rmse_tol)
    else:  # use absolute tolerance
        self.target_rmse_tol = float(abs_tol) / norm.ppf(1 - alpha / 2)
    self.rmse_tol = self.target_rmse_tol  # user-facing attribute; never mutated after __init__
    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._active_trace = None
    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"

integrate

integrate(resume=None)

Run (or continue) the continuation-MLQMC integration.

Parameters:

Name	Type	Description	Default
`resume`	`Data`	Checkpoint returned by a previous `integrate()` call. The new tolerance may be tighter or looser than the one used when the checkpoint was created. With a tighter tolerance the algorithm picks up the tolerance ladder from `max(checkpoint_rmse_tol, target_rmse_tol)` and continues down to `target_rmse_tol`. With a looser tolerance the first step immediately converges on the existing samples and no additional ladder steps are needed.	`None`

Returns:

Name	Type	Description
`tuple`	`tuple`	`(solution, data)`.

Source code in qmcpy/stopping_criterion/cub_mlqmc_cont.py

def integrate(self, resume=None) -> tuple:
    """Run (or continue) the continuation-MLQMC integration.

    Args:
        resume (Data, optional): Checkpoint returned by a previous
            ``integrate()`` call.  The new tolerance may be tighter *or*
            looser than the one used when the checkpoint was created.
            With a tighter tolerance the algorithm picks up the tolerance
            ladder from ``max(checkpoint_rmse_tol, target_rmse_tol)`` and
            continues down to ``target_rmse_tol``.  With a looser tolerance
            the first step immediately converges on the existing samples
            and no additional ladder steps are needed.

    Returns:
        tuple: ``(solution, data)``.
    """
    self._active_t_start = t_start = time()
    self._active_trace = trace = self._make_trace_logger()
    self._active_resume_provenance = resume_provenance = self._capture_resume_provenance(resume)
    try:
        data = self._prepare_resume_data(resume, self._validate_resume, self._restore_resume_state)
        if data is not None:
            step_tol = max(getattr(data, 'rmse_tol', self.target_rmse_tol), self.target_rmse_tol)
            data.rmse_tol = step_tol
            self._set_elapsed_time(data, 0.0, resume_provenance=resume_provenance)
            trace.resume(data, step_value=int(data.levels))
            if step_tol <= self.target_rmse_tol:
                # Same or looser tolerance: check convergence at target_tol
                self._integrate(
                    data,
                    skip_level_reset=True,
                    step_tol=step_tol,
                )
            else:
                # Tighter tolerance: skip to ladder, preserving level structure for first step
                first = True
                for t in range(self.n_tols):
                    next_tol = self.inflate ** (self.n_tols - t - 1) * self.target_rmse_tol
                    if next_tol < step_tol:
                        self._integrate(
                            data,
                            skip_level_reset=first,
                            step_tol=next_tol,
                        )
                        first = False
        else:
            data = self._construct_data()
            # Loop over coarser tolerances
            for t in range(self.n_tols):
                step_tol = self.inflate ** (self.n_tols - t - 1) * self.target_rmse_tol
                self._integrate(
                    data,
                    step_tol=step_tol,
                )
        self._finalize_integration_data(
            data, time() - t_start, resume_provenance=resume_provenance
        )
        trace.finalize()
        data.iteration_history = getattr(self, "iteration_history", None)
        data.history_df = getattr(self, "history_df", None)
        return data.solution, data
    finally:
        self._active_trace = None
        self._active_t_start = None
        self._active_resume_provenance = None

`CubMLQMC`

Bases: AbstractCubMLQMC

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

def __init__(
    self,
    integrand,
    abs_tol=0.05,
    rmse_tol=None,
    n_init=256,
    n_limit=1e10,
    alpha=0.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"

integrate

integrate(resume=None)

Run (or continue) the MLQMC integration.

Parameters:

Name	Type	Description	Default
`resume`	`Data`	Checkpoint returned by a previous `integrate()` call. The new tolerance may be tighter or looser than the one used when the checkpoint was created. With a tighter tolerance the algorithm draws additional samples from where it left off. With a looser tolerance the existing samples already satisfy the requirement and the method returns immediately with no new sampling.	`None`

Returns:

Name	Type	Description
`tuple`	`tuple`	`(solution, data)`.

Source code in qmcpy/stopping_criterion/cub_mlqmc.py

def integrate(self, resume=None) -> tuple:
    """Run (or continue) the MLQMC integration.

    Args:
        resume (Data, optional): Checkpoint returned by a previous
            ``integrate()`` call.  The new tolerance may be tighter *or*
            looser than the one used when the checkpoint was created.
            With a tighter tolerance the algorithm draws additional samples
            from where it left off.  With a looser tolerance the existing
            samples already satisfy the requirement and the method returns
            immediately with no new sampling.

    Returns:
        tuple: ``(solution, data)``.
    """
    t_start = time()
    resume_provenance = self._capture_resume_provenance(resume)
    trace = self._make_trace_logger()
    data = self._prepare_resume_data(resume, self._validate_resume, self._restore_resume_state)
    if data is not None:
        data.rmse_estimate = np.sqrt(data.var_level[:data.levels].sum() + data.bias_estimate**2)
        self._set_elapsed_time(data, 0.0, resume_provenance=resume_provenance)
        trace.resume(data, step_value=int(data.levels))
    if data is None:
        data = self._construct_data()
    self._run_integrate_loop(
        data,
        self.update_data,
        trace=trace,
        t_start=t_start,
        resume_provenance=resume_provenance,
    )
    self._finalize_integration_data(
        data, time() - t_start, resume_provenance=resume_provenance
    )
    trace.finalize()
    data.iteration_history = getattr(self, "iteration_history", None)
    data.history_df = getattr(self, "history_df", None)
    return data.solution, data

Multilevel IID MC

`CubMLMCCont`

Bases: AbstractCubMLMC

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

def __init__(
    self,
    integrand,
    abs_tol=0.05,
    rmse_tol=None,
    n_init=256,
    n_limit=1e10,
    inflate=100 ** (1 / 9),
    alpha=0.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_rmse_tol = float(rmse_tol)
    else:  # use absolute tolerance
        self.target_rmse_tol = float(abs_tol) / norm.ppf(1 - alpha / 2)
    self.rmse_tol = self.target_rmse_tol  # user-facing attribute; never mutated after __init__
    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._active_trace = None
    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,
    )

integrate

integrate(resume=None)

Run (or continue) the continuation-MLMC integration.

Parameters:

Name	Type	Description	Default
`resume`	`Data`	Checkpoint returned by a previous `integrate()` call. The new tolerance may be tighter or looser than the one used when the checkpoint was created. With a tighter tolerance the algorithm picks up the tolerance ladder from `max(checkpoint_rmse_tol, target_rmse_tol)` and continues down to `target_rmse_tol`. With a looser tolerance the first step immediately converges on the existing samples and no additional ladder steps are needed.	`None`

Returns:

Name	Type	Description
`tuple`	`tuple`	`(solution, data)`.

Source code in qmcpy/stopping_criterion/cub_mlmc_cont.py

def integrate(self, resume=None) -> tuple:
    """Run (or continue) the continuation-MLMC integration.

    Args:
        resume (Data, optional): Checkpoint returned by a previous
            ``integrate()`` call.  The new tolerance may be tighter *or*
            looser than the one used when the checkpoint was created.
            With a tighter tolerance the algorithm picks up the tolerance
            ladder from ``max(checkpoint_rmse_tol, target_rmse_tol)`` and
            continues down to ``target_rmse_tol``.  With a looser tolerance
            the first step immediately converges on the existing samples
            and no additional ladder steps are needed.

    Returns:
        tuple: ``(solution, data)``.
    """
    self._active_t_start = t_start = time()
    self._active_trace = trace = self._make_trace_logger()
    self._active_resume_provenance = resume_provenance = self._capture_resume_provenance(resume)
    try:
        data = self._prepare_resume_data(resume, self._validate_resume, self._restore_resume_state)
        replay_snapshots = None
        replay_iter_count = None
        if self._can_replay_resume_exactly(data):
            replay_data, replay_snapshots, replay_iter_count = self._replay_resume_exactly(
                data, t_start=t_start, resume_provenance=resume_provenance
            )
            if replay_data is not None:
                data = replay_data
        if resume is not None:
            checkpoint = self._prepare_resume_data(
                resume, self._validate_resume, self._restore_resume_state
            )
            step_tol = max(
                getattr(checkpoint, "rmse_tol", self.target_rmse_tol),
                self.target_rmse_tol,
            )
            checkpoint.rmse_tol = step_tol
            if replay_iter_count is not None:
                checkpoint._iter_count = replay_iter_count
            self._set_elapsed_time(checkpoint, 0.0, resume_provenance=resume_provenance)
            trace.resume(checkpoint, step_value=int(checkpoint.levels + 1))
        if replay_snapshots is not None:
            for snapshot in replay_snapshots:
                trace.iteration(snapshot, step_value=int(snapshot.levels + 1))
            if replay_snapshots and trace.enabled:
                data._iter_count = trace.iter_count
        elif data is not None:
            step_tol = max(getattr(data, 'rmse_tol', self.target_rmse_tol), self.target_rmse_tol)
            data.rmse_tol = step_tol
            if step_tol <= self.target_rmse_tol:
                # Same or looser tolerance: check convergence at target_tol
                self._integrate(
                    data,
                    skip_level_reset=True,
                    step_tol=step_tol,
                )
            else:
                # Tighter tolerance: skip to ladder, preserving level structure for first step
                first = True
                for t in range(self.n_tols):
                    next_tol = self.inflate ** (self.n_tols - t - 1) * self.target_rmse_tol
                    if next_tol < step_tol:
                        self._integrate(
                            data,
                            skip_level_reset=first,
                            step_tol=next_tol,
                        )
                        first = False
        else:
            data = self._construct_data()
            # Loop over coarser tolerances
            for t in range(self.n_tols):
                step_tol = self.inflate ** (self.n_tols - t - 1) * self.target_rmse_tol
                self._integrate(
                    data,
                    step_tol=step_tol,
                )
        self._finalize_integration_data(
            data, time() - t_start, resume_provenance=resume_provenance
        )
        trace.finalize()
        data.iteration_history = getattr(self, "iteration_history", None)
        data.history_df = getattr(self, "history_df", None)
        return data.solution, data
    finally:
        self._active_trace = None
        self._active_t_start = None
        self._active_resume_provenance = None

`CubMLMC`

Bases: AbstractCubMLMC

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

def __init__(
    self,
    integrand,
    abs_tol=0.05,
    rmse_tol=None,
    n_init=256,
    n_limit=1e10,
    alpha=0.01,
    levels_min=2,
    levels_max=10,
    alpha0=-1.0,
    beta0=-1.0,
    gamma0=-1.0,
):
    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,
    )

integrate

integrate(resume=None)

Run (or continue) the MLMC integration.

Parameters:

Name	Type	Description	Default
`resume`	`Data`	Checkpoint returned by a previous `integrate()` call. The new tolerance may be tighter or looser than the one used when the checkpoint was created. With a tighter tolerance the algorithm draws additional samples from where it left off. With a looser tolerance the existing samples already satisfy the requirement and the method returns immediately with no new sampling.	`None`

Returns:

Name	Type	Description
`tuple`	`tuple`	`(solution, data)`.

Source code in qmcpy/stopping_criterion/cub_mlmc.py

def integrate(self, resume=None) -> tuple:
    """Run (or continue) the MLMC integration.

    Args:
        resume (Data, optional): Checkpoint returned by a previous
            ``integrate()`` call.  The new tolerance may be tighter *or*
            looser than the one used when the checkpoint was created.
            With a tighter tolerance the algorithm draws additional samples
            from where it left off.  With a looser tolerance the existing
            samples already satisfy the requirement and the method returns
            immediately with no new sampling.

    Returns:
        tuple: ``(solution, data)``.
    """
    t_start = time()
    resume_provenance = self._capture_resume_provenance(resume)
    trace = self._make_trace_logger()
    data = self._prepare_resume_data(resume, self._validate_resume, self._restore_resume_state)
    replay_snapshots = None
    replay_iter_count = None
    if self._can_replay_resume_exactly(data):
        replay_data, replay_snapshots, replay_iter_count = self._replay_resume_exactly(
            data, t_start=t_start, resume_provenance=resume_provenance
        )
        if replay_data is not None:
            data = replay_data
    if resume is not None:
        checkpoint = self._prepare_resume_data(
            resume, self._validate_resume, self._restore_resume_state
        )
        checkpoint.rmse_estimate = np.sqrt(
            (checkpoint.var_level / checkpoint.n_level[: len(checkpoint.var_level)]).sum()
        )
        if replay_iter_count is not None:
            checkpoint._iter_count = replay_iter_count
        self._set_elapsed_time(checkpoint, 0.0, resume_provenance=resume_provenance)
        trace.resume(checkpoint, step_value=int(checkpoint.levels + 1))
    if replay_snapshots is not None:
        for snapshot in replay_snapshots:
            trace.iteration(snapshot, step_value=int(snapshot.levels + 1))
        if replay_snapshots and trace.enabled:
            data._iter_count = trace.iter_count
    else:
        if data is None:
            data = self._construct_data()
        self._run_integrate_loop(
            data,
            self._update_data,
            trace=trace,
            t_start=t_start,
            resume_provenance=resume_provenance,
        )
    data.rmse_tol = self.rmse_tol
    data.levels += 1
    self._finalize_integration_data(
        data, time() - t_start, resume_provenance=resume_provenance
    )
    trace.finalize()
    data.iteration_history = getattr(self, "iteration_history", None)
    data.history_df = getattr(self, "history_df", None)
    return data.solution, data

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:

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

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

Stopping Criteria

UML Overview

AbstractStoppingCriterion

integrate

get_iteration_log

format_iteration_log

print_iteration_log

set_tolerance

QMC

CubQMCNetG

CubQMCLatticeG

CubQMCBayesNetG

CubQMCBayesLatticeG

CubQMCRepStudentT

IID MC

CubMCCLTVec

CubMCG

CubMCCLT

Multilevel QMC

CubMLQMCCont

integrate

CubMLQMC

integrate

Multilevel IID MC

CubMLMCCont

integrate

CubMLMC

integrate

Failure Probability

PFGPCI

UML Specific

`AbstractStoppingCriterion`

`CubQMCNetG`

`CubQMCLatticeG`

`CubQMCBayesNetG`

`CubQMCBayesLatticeG`

`CubQMCRepStudentT`

`CubMCCLTVec`

`CubMCG`

`CubMCCLT`

`CubMLQMCCont`

`CubMLQMC`

`CubMLMCCont`

`CubMLMC`

`PFGPCI`