Integrands

UML Overview

AbstractIntegrand

Bases: object

Parameters:

Name	Type	Description	Default
dimension_indv	tuple	Individual solution shape.	required
dimension_comb	tuple	Combined solution shape.	required
parallel	int	Parallelization flag. When parallel = 0 or parallel = 1 then function evaluation is done in serial fashion. parallel > 1 specifies the number of processes used by multiprocessing.Pool or multiprocessing.pool.ThreadPool. Setting parallel=True is equivalent to parallel = os.cpu_count().	required
threadpool	bool	When parallel > 1: Setting threadpool = True will use multiprocessing.pool.ThreadPool. Setting threadpool = False will use setting multiprocessing.Pool.	False

Source code in qmcpy/integrand/abstract_integrand.py

def __init__(self, dimension_indv, dimension_comb, parallel, threadpool=False):
    r"""
    Args:
        dimension_indv (tuple): Individual solution shape.
        dimension_comb (tuple): Combined solution shape. 
        parallel (int): Parallelization flag. 

            - When `parallel = 0` or `parallel = 1` then function evaluation is done in serial fashion.
            - `parallel > 1` specifies the number of processes used by `multiprocessing.Pool` or `multiprocessing.pool.ThreadPool`.

            Setting `parallel=True` is equivalent to `parallel = os.cpu_count()`.
        threadpool (bool): When `parallel > 1`: 

            - Setting `threadpool = True` will use `multiprocessing.pool.ThreadPool`.
            - Setting `threadpool = False` will use `setting multiprocessing.Pool`.
    """
    prefix = 'A concrete implementation of AbstractIntegrand must have '
    self.d = self.true_measure.d
    self.d_indv = (dimension_indv,) if isinstance(dimension_indv,int) else tuple(dimension_indv)
    self.d_comb = (dimension_comb,) if isinstance(dimension_comb,int) else tuple(dimension_comb)
    cpus = os.cpu_count()
    self.parallel = cpus if parallel is True else int(parallel)
    self.parallel = 0 if self.parallel==1 else self.parallel
    self.threadpool = threadpool
    if not (hasattr(self, 'sampler') and isinstance(self.sampler,(AbstractTrueMeasure,AbstractDiscreteDistribution))):
        raise ParameterError(prefix + 'self.sampler, a AbstractTrueMeasure or AbstractDiscreteDistribution instance')
    if not (hasattr(self, 'true_measure') and isinstance(self.true_measure,AbstractTrueMeasure)):
        raise ParameterError(prefix + 'self.true_measure, a AbstractTrueMeasure instance')
    if not hasattr(self,'parameters'):
        self.parameters = []
    if not hasattr(self,'multilevel'):
        self.multilevel = False
    assert isinstance(self.multilevel,bool)
    if not hasattr(self,'max_level'):
        self.max_level = np.inf
    if not hasattr(self,'discrete_distrib'):
        self.discrete_distrib = self.true_measure.discrete_distrib
    if self.true_measure.transform!=self.true_measure and \
       not (self.true_measure.range==self.true_measure.transform.range).all():
        raise ParameterError("The range of the composed transform is not compatible with this true measure")
    self.EPS = np.finfo(float).eps

__call__

__call__(n=None, n_min=None, n_max=None, warn=True)

• If just n is supplied, generate samples from the sequence at indices 0,...,n-1.
• If n_min and n_max are supplied, generate samples from the sequence at indices n_min,...,n_max-1.
• If n and n_min are supplied, then generate samples from the sequence at indices n,...,n_min-1.

Parameters:

Name	Type	Description	Default
n	Union[None, int]	Number of points to generate.	None
n_min	Union[None, int]	Starting index of sequence.	None
n_max	Union[None, int]	Final index of sequence.	None
warn	bool	If False, disable warnings when generating samples.	True

Returns:

Name	Type	Description
t	ndarray	Samples from the sequence. If replications is None then this will be of size (n_max-n_min) \(\times\) dimension If replications is a positive int, then t will be of size replications \(\times\) (n_max-n_min) \(\times\) dimension
weights	ndarray	Only returned when return_weights=True. The Jacobian weights for the transformation

Source code in qmcpy/integrand/abstract_integrand.py

def __call__(self, n=None, n_min=None, n_max=None, warn=True):
    r"""
    - If just `n` is supplied, generate samples from the sequence at indices 0,...,`n`-1.
    - If `n_min` and `n_max` are supplied, generate samples from the sequence at indices `n_min`,...,`n_max`-1.
    - If `n` and `n_min` are supplied, then generate samples from the sequence at indices `n`,...,`n_min`-1.

    Args:
        n (Union[None,int]): Number of points to generate.
        n_min (Union[None,int]): Starting index of sequence.
        n_max (Union[None,int]): Final index of sequence.
        warn (bool): If `False`, disable warnings when generating samples.

    Returns:
        t (np.ndarray): Samples from the sequence. 

            - If `replications` is `None` then this will be of size (`n_max`-`n_min`) $\times$ `dimension` 
            - If `replications` is a positive int, then `t` will be of size `replications` $\times$ (`n_max`-`n_min`) $\times$ `dimension` 
        weights (np.ndarray): Only returned when `return_weights=True`. The Jacobian weights for the transformation
    """
    return self.gen_samples(n=n,n_min=n_min,n_max=n_max,warn=warn)

g

g(t, *args, **kwargs)

Abstract method implementing the integrand as a function of the true measure.

Parameters:

Name	Type	Description	Default
t	ndarray	Inputs with shape (*batch_shape, d).	required
args	tuple	positional arguments to g.	()
kwargs	dict	keyword arguments to g. Some algorithms will additionally try to pass in a compute_flags keyword argument. This np.ndarray are flags indicating which outputs require evaluation. For example, if the vector function has 3 outputs and compute_flags = [False, True, False], then the function is only required to evaluate the second output and may leave the remaining outputs as np.nan values, i.e., the outputs corresponding to compute_flags which are False will not be used in the computation.	{}

Returns:

Name	Type	Description
y	ndarray	function evaluations with shape (*batch_shape, *dimension_indv) where dimension_indv is the shape of the function outputs.

Source code in qmcpy/integrand/abstract_integrand.py

def g(self, t, *args, **kwargs):
    r"""
    *Abstract method* implementing the integrand as a function of the true measure.

    Args:
        t (np.ndarray): Inputs with shape `(*batch_shape, d)`.
        args (tuple): positional arguments to `g`.
        kwargs (dict): keyword arguments to `g`. 

            Some algorithms will additionally try to pass in a `compute_flags` keyword argument. 
            This `np.ndarray` are flags indicating which outputs require evaluation.  
            For example, if the vector function has 3 outputs and `compute_flags = [False, True, False]`, 
            then the function is only required to evaluate the second output and may leave the remaining outputs as `np.nan` values, 
            i.e., the outputs corresponding to `compute_flags` which are `False` will not be used in the computation.

    Returns:
        y (np.ndarray): function evaluations with shape `(*batch_shape, *dimension_indv)` where `dimension_indv` is the shape of the function outputs. 
    """
    raise MethodImplementationError(self, 'g')

f

f(x, *args, **kwargs)

Function to evaluate the transformed integrand as a function of the discrete distribution.
Automatically applies the transformation determined by the true measure.

Parameters:

Name	Type	Description	Default
x	ndarray	Inputs with shape (*batch_shape, d).	required
args	tuple	positional arguments to g.	()
kwargs	dict	keyword arguments to g. Some algorithms will additionally try to pass in a compute_flags keyword argument. This np.ndarray are flags indicating which outputs require evaluation. For example, if the vector function has 3 outputs and compute_flags = [False, True, False], then the function is only required to evaluate the second output and may leave the remaining outputs as np.nan values, i.e., the outputs corresponding to compute_flags which are False will not be used in the computation. The keyword argument periodization_transform, a string, specifies a periodization transform. Options are: False: No periodizing transform, \(\psi(x) = x\). 'BAKER': Baker tansform \(\psi(x) = 1-2\lvert x-1/2 \rvert\). 'C0': \(C^0\) transform \(\psi(x) = 3x^2-2x^3\). 'C1': \(C^1\) transform \(\psi(x) = x^3(10-15x+6x^2)\). 'C1SIN': Sidi \(C^1\) transform \(\psi(x) = x-\sin(2 \pi x)/(2 \pi)\). 'C2SIN': Sidi \(C^2\) transform \(\psi(x) = (8-9 \cos(\pi x)+\cos(3 \pi x))/16\). 'C3SIN': Sidi \(C^3\) transform \(\psi(x) = (12\pi x-8\sin(2 \pi x) + \sin(4 \pi x))/(12 \pi)\).	{}

Returns:

Name	Type	Description
y	ndarray	function evaluations with shape (*batch_shape, *dimension_indv) where dimension_indv is the shape of the function outputs.

Source code in qmcpy/integrand/abstract_integrand.py

def f(self, x, *args, **kwargs):
    r"""
    Function to evaluate the transformed integrand as a function of the discrete distribution.  
    Automatically applies the transformation determined by the true measure. 

    Args:
        x (np.ndarray): Inputs with shape `(*batch_shape, d)`.
        args (tuple): positional arguments to `g`.
        kwargs (dict): keyword arguments to `g`. 

            Some algorithms will additionally try to pass in a `compute_flags` keyword argument. 
            This `np.ndarray` are flags indicating which outputs require evaluation.  
            For example, if the vector function has 3 outputs and `compute_flags = [False, True, False]`, 
            then the function is only required to evaluate the second output and may leave the remaining outputs as `np.nan` values, 
            i.e., the outputs corresponding to `compute_flags` which are `False` will not be used in the computation.

            The keyword argument `periodization_transform`, a string, specifies a periodization transform. 
            Options are: 

            - `False`: No periodizing transform, $\psi(x) = x$. 
            - `'BAKER'`: Baker tansform $\psi(x) = 1-2\lvert x-1/2 \rvert$.
            - `'C0'`: $C^0$ transform $\psi(x) = 3x^2-2x^3$.
            - `'C1'`: $C^1$ transform $\psi(x) = x^3(10-15x+6x^2)$.
            - `'C1SIN'`: Sidi $C^1$ transform $\psi(x) = x-\sin(2 \pi x)/(2 \pi)$. 
            - `'C2SIN'`: Sidi $C^2$ transform $\psi(x) = (8-9 \cos(\pi x)+\cos(3 \pi x))/16$.
            - `'C3SIN'`: Sidi $C^3$ transform $\psi(x) = (12\pi x-8\sin(2 \pi x) + \sin(4 \pi x))/(12 \pi)$.

    Returns:
        y (np.ndarray): function evaluations with shape `(*batch_shape, *dimension_indv)` where `dimension_indv` is the shape of the function outputs. 
    """
    if "periodization_transform" in kwargs:
        periodization_transform = kwargs["periodization_transform"]
        del kwargs["periodization_transform"]
    else:
        periodization_transform = "None"
    periodization_transform = str(periodization_transform).upper()
    batch_shape = tuple(x.shape[:-1])
    d_indv_ndim = len(self.d_indv)
    if periodization_transform=="NONE": periodization_transform = "FALSE"
    if self.discrete_distrib.mimics != 'StdUniform' and periodization_transform!='NONE':
        raise ParameterError('''
            Applying a periodization transform currently requires a discrete distribution 
            that mimics a standard uniform measure.''')
    if periodization_transform == 'FALSE':
        xp = x
        wp = np.ones(batch_shape,dtype=float)
    elif periodization_transform == 'BAKER':
        xp = 1-2*abs(x-1/2)
        wp = np.ones(batch_shape,dtype=float)
    elif periodization_transform == 'C0':
        xp = 3*x**2-2*x**3
        wp = np.prod(6*x*(1-x),-1)
    elif periodization_transform == 'C1':
        xp = x**3*(10-15*x+6*x**2)
        wp = np.prod(30*x**2*(1-x)**2,-1)
    elif periodization_transform == 'C1SIN':
        xp = x - np.sin(2*np.pi*x)/(2*np.pi)
        wp = np.prod(2*np.sin(np.pi*x)**2,-1)
    elif periodization_transform == 'C2SIN':
        xp = (8-9*np.cos(np.pi*x)+np.cos(3*np.pi*x))/16
        wp = np.prod((9*np.sin(np.pi*x)*np.pi-np.sin(3*np.pi*x)*3*np.pi)/16,-1)
    elif periodization_transform=='C3SIN':
        xp = (12*np.pi*x-8*np.sin(2*np.pi*x)+np.sin(4*np.pi*x))/(12*np.pi)
        wp = np.prod((12*np.pi-8*np.cos(2*np.pi*x)*2*np.pi+np.sin(4*np.pi*x)*4*np.pi)/(12*np.pi),-1)
    else:
        raise ParameterError("The %s periodization transform is not implemented"%periodization_transform)
    if periodization_transform in ['C1','C1SIN','C2SIN','C3SIN']:
        xp[xp<=0] = self.EPS
        xp[xp>=1] = 1-self.EPS
    assert wp.shape==batch_shape
    assert xp.shape==x.shape
    # function evaluation with chain rule
    i = (None,)*d_indv_ndim+(...,)
    if self.true_measure==self.true_measure.transform:
        # jacobian*weight/pdf will cancel so f(x) = g(\Psi(x))
        xtf = self.true_measure._jacobian_transform_r(xp,return_weights=False) # get transformed samples, equivalent to self.true_measure._transform_r(x)
        assert xtf.shape==xp.shape
        y = self._g(xtf,*args,**kwargs)
    else: # using importance sampling --> need to compute pdf, jacobian(s), and weight explicitly
        pdf = self.discrete_distrib.pdf(xp) # pdf of samples
        assert pdf.shape==batch_shape
        xtf,jacobians = self.true_measure.transform._jacobian_transform_r(xp,return_weights=True) # compute recursive transform+jacobian
        assert xtf.shape==xp.shape 
        assert jacobians.shape==batch_shape
        weight = self.true_measure._weight(xtf) # weight based on the true measure
        assert weight.shape==batch_shape
        gvals = self._g(xtf,*args,**kwargs)
        assert gvals.shape==(self.d_indv+batch_shape)
        y = gvals*weight[i]/pdf[i]*jacobians[i]
    assert y.shape==(self.d_indv+batch_shape)
    # account for periodization weight
    y = y*wp[i]
    assert y.shape==(self.d_indv+batch_shape)
    return y

bound_fun

bound_fun(bound_low, bound_high)

Compute the bounds on the combined function based on bounds for the individual functions.

Defaults to the identity where we essentially do not combine integrands, but instead integrate each function individually.

Parameters:

Name	Type	Description	Default
bound_low	ndarray	Lower bounds on individual estimates with shape integrand.d_indv.	required
bound_high	ndarray	Upper bounds on individual estimates with shape integrand.d_indv.	required

Returns:

Name	Type	Description
comb_bound_low	ndarray	Lower bounds on combined estimates with shape integrand.d_comb.
comb_bound_high	ndarray	Upper bounds on combined estimates with shape integrand.d_comb.

Source code in qmcpy/integrand/abstract_integrand.py

def bound_fun(self, bound_low, bound_high):
    """
    Compute the bounds on the combined function based on bounds for the
    individual functions.  

    Defaults to the identity where we essentially
    do not combine integrands, but instead integrate each function
    individually.

    Args:
        bound_low (np.ndarray): Lower bounds on individual estimates with shape `integrand.d_indv`.
        bound_high (np.ndarray): Upper bounds on individual estimates with shape `integrand.d_indv`.

    Returns:
        comb_bound_low (np.ndarray): Lower bounds on combined estimates with shape `integrand.d_comb`.
        comb_bound_high (np.ndarray): Upper bounds on combined estimates with shape `integrand.d_comb`.
    """
    if self.d_indv!=self.d_comb:
        raise ParameterError('''
            Set bound_fun explicitly. 
            The default bound_fun is the identity map. 
            Since the individual solution dimensions d_indv = %s does not equal the combined solution dimensions d_comb = %d, 
            QMCPy cannot infer a reasonable bound function.'''%(str(self.d_indv),str(self.d_comb)))
    return bound_low,bound_high

dependency

dependency(comb_flags)

Takes a vector of indicators of weather of not the error bound is satisfied for combined integrands and returns flags for individual integrands.

For example, if we are taking the ratio of 2 individual integrands, then getting comb_flags=True means the ratio has not been approximated to within the tolerance, so the dependency function should return indv_flags=[True,True] indicating that both the numerator and denominator integrands need to be better approximated.

Parameters:

Name	Type	Description	Default
comb_flags	ndarray	Flags of shape integrand.d_comb indicating whether the combined outputs are insufficiently approximated.	required

Returns:

Name	Type	Description
indv_flags	ndarray	Flags of shape integrand.d_indv indicating whether the individual integrands require additional sampling.

Source code in qmcpy/integrand/abstract_integrand.py

def dependency(self, comb_flags):
    """
    Takes a vector of indicators of weather of not the error bound is satisfied for combined integrands and returns flags for individual integrands.  

    For example, if we are taking the ratio of 2 individual integrands, then getting `comb_flags=True` means the ratio
    has not been approximated to within the tolerance, so the dependency function should return `indv_flags=[True,True]`
    indicating that both the numerator and denominator integrands need to be better approximated.

    Args:
        comb_flags (np.ndarray): Flags of shape `integrand.d_comb` indicating whether the combined outputs are insufficiently approximated.

    Returns:
        indv_flags (np.ndarray): Flags of shape `integrand.d_indv` indicating whether the individual integrands require additional sampling. 
    """
    return comb_flags if self.d_indv==self.d_comb else np.tile((comb_flags==False).any(),self.d_indv)

spawn

spawn(levels)

Spawn new instances of the current integrand at different levels with new seeds. Used by multi-level QMC algorithms which require integrands at multiple levels.

Note

Use replications instead of using spawn when possible, e.g., when spawning copies which all have the same level.

Parameters:

Name	Type	Description	Default
levels	ndarray	Levels at which to spawn new integrands.	required

Returns:

Name	Type	Description
spawned_integrand	list	Integrands with new true measures and discrete distributions.

Source code in qmcpy/integrand/abstract_integrand.py

def spawn(self, levels):
    r"""
    Spawn new instances of the current integrand at different levels with new seeds.
    Used by multi-level QMC algorithms which require integrands at multiple levels.

    Note:
        Use `replications` instead of using `spawn` when possible, e.g., when spawning copies which all have the same level.

    Args:
        levels (np.ndarray): Levels at which to spawn new integrands.

    Returns:
        spawned_integrand (list): Integrands with new true measures and discrete distributions.
    """
    levels = np.array([levels]) if np.isscalar(levels) else np.array(levels)
    if (levels>self.max_level).any():
        raise ParameterError("requested spawn level exceeds max level")
    n_levels = len(levels)
    new_dims = np.array([self.dimension_at_level(level) for level in levels])
    tm_spawns = self.sampler.spawn(s=n_levels,dimensions=new_dims) 
    spawned_integrand = [None]*n_levels
    for l,level in enumerate(levels):
        spawned_integrand[l] = self._spawn(level,tm_spawns[l])
    return spawned_integrand

dimension_at_level

dimension_at_level(level)

Abstract method which returns the dimension of the generator required for a given level.

Note

Only used for multilevel problems.

Parameters:

Name	Type	Description	Default
level	int	Level at which to return the dimension.	required

Returns:

Name	Type	Description
d	int	Dimension at the given input level.

Source code in qmcpy/integrand/abstract_integrand.py

def dimension_at_level(self, level):
    """
    *Abstract method* which returns the dimension of the generator required for a given level.

    Note:
        Only used for multilevel problems.

    Args:
        level (int): Level at which to return the dimension.

    Returns:
        d (int): Dimension at the given input level. 
    """
    return self.d

FinancialOption

Bases: AbstractIntegrand

Financial options.

• Start price \(S_0\)
• Strike price \(K\)
• Interest rate \(r\)
• Volatility \(\sigma\)
• Equidistant monitoring times \(\boldsymbol{\tau} = (\tau_1,\dots,\tau_d)^T\) with \(\tau_d\) the final (exercise) time and \(\tau_j = \tau_d j/d\).

Define the geometric brownian motion as

\[\boldsymbol{S}(\boldsymbol{t}) = S_0 e^{(r-\sigma^2/2)\boldsymbol{\tau}+\sigma\boldsymbol{t}}, \qquad \boldsymbol{T} \sim \mathcal{N}(\boldsymbol{0},\mathsf{\Sigma})\]

where \(\boldsymbol{T}\) is a standard Brownian motion so \(\mathsf{\Sigma} = \left(\min\{\tau_j,\tau_{j'}\}\right)_{j,j'=1}^d\).

The discounted payoff is

\[g(\boldsymbol{t}) = P(\boldsymbol{S}(\boldsymbol{t}))e^{-r \tau_d}\]

where the payoff function \(P\) will be defined depending on the option.

Below we wil luse \(S_{-1}\) to denote the final element of \(\boldsymbol{S}\), the value of the path at exercise time.

European Options

European Call and Put Options have respective payoffs

\[P(\boldsymbol{S}) = \max\{S_{-1}-K,0\}, \qquad P(\boldsymbol{S}) = \max\{K-S_{-1},0\}.\]

Asian Options

An asian option considers the average value of an asset path across time. We use the trapezoidal rule to approximate either the arithmetic mean by

\[A(\boldsymbol{S}) = \frac{1}{d}\left[\frac{1}{2} S_0 + \sum_{j=1}^{d-1} S_j + \frac{1}{2} S_{-1}\right]\]

or the geometric mean by

\[A(\boldsymbol{S}) = \left[\sqrt{S_0} \prod_{j=1}^{d-1} S_j \sqrt{S_{-1}}\right]^{1/d}.\]

Asian Call and Put Option have respective payoffs

\[P(\boldsymbol{S}) = \max\{A(\boldsymbol{S})-K,0\}, \qquad P(\boldsymbol{S}) = \max\{K-A(\boldsymbol{S}),0\}.\]

Barrier Options

• Barrier \(B\).

In options are activate when the path crosses the barrier \(B\), while out options are activated only if the path never crosses the barrier \(B\). An up option satisfies \(S_0<B\) while a down option satisfies \(S_0>B\), both indicating the direction of the barrier from the start price.

Barrier Up-In Call and Put Options have respective payoffs

\[P(\boldsymbol{S}) = \begin{cases} \max\{S_{-1})-K,0\}, & \text{any } \boldsymbol{S} \geq B \\ 0, & \mathrm{otherwise} \end{cases}, \qquad P(\boldsymbol{S}) = \begin{cases} \max\{K-S_{-1}),0\}, & \text{any } \boldsymbol{S} \geq B \\ 0, & \mathrm{otherwise} \end{cases}.\]

Barrier Up-Out Call and Put Options have respective payoffs

\[P(\boldsymbol{S}) = \begin{cases} \max\{S_{-1})-K,0\}, & \text{all } \boldsymbol{S} < B \\ 0, & \mathrm{otherwise} \end{cases}, \qquad P(\boldsymbol{S}) = \begin{cases} \max\{K-S_{-1}),0\}, & \text{all } \boldsymbol{S} < B \\ 0, & \mathrm{otherwise} \end{cases}.\]

Barrier Down-In Call and Put Options have respective payoffs

\[P(\boldsymbol{S}) = \begin{cases} \max\{S_{-1})-K,0\}, & \text{any } \boldsymbol{S} \leq B \\ 0, & \mathrm{otherwise} \end{cases}, \qquad P(\boldsymbol{S}) = \begin{cases} \max\{K-S_{-1}),0\}, & \text{any } \boldsymbol{S} \leq B \\ 0, & \mathrm{otherwise} \end{cases}.\]

Barrier Down-Out Call and Put Options have respective payoffs

\[P(\boldsymbol{S}) = \begin{cases} \max\{S_{-1})-K,0\}, & \text{all } \boldsymbol{S} > B \\ 0, & \mathrm{otherwise} \end{cases}, \qquad P(\boldsymbol{S}) = \begin{cases} \max\{K-S_{-1}),0\}, & \text{all } \boldsymbol{S} > B \\ 0, & \mathrm{otherwise} \end{cases}.\]

Lookback Options

Lookback Call and Put Options have respective payoffs

\[P(\boldsymbol{S}) = S_{-1}-\min(\boldsymbol{S}), \qquad P(\boldsymbol{S}) = \max(\boldsymbol{S})-S_{-1}.\]

Digital Option

• Payout \(\rho\).

Digital Call and Put Options have respective payoffs

\[P(\boldsymbol{S}) = \begin{cases} \rho, & S_{-1} \geq K \\ 0, & \mathrm{otherwise} \end{cases}, \qquad P(\boldsymbol{S}) = \begin{cases} \rho, & S_{-1} \leq K \\ 0, & \mathrm{otherwise} \end{cases}.\]

Multilevel Options

• Initial level \(\ell_0 \geq 0\).
• Level \(\ell \geq \ell_0\).

Let \(\boldsymbol{S}_\mathrm{fine}=\boldsymbol{S}\) be the fine full path. For \(\ell>\ell_0\) write the coarse path as \(\boldsymbol{S}_\mathrm{coarse} = (S_j)_{j \text{ even}}\) which only considers every other element of \(\boldsymbol{S}\). In this multilevel setting the payoff is

\[P_\ell(\boldsymbol{S}) = \begin{cases} P(\boldsymbol{S}_\mathrm{fine}), & \ell = \ell_0, \\ P(\boldsymbol{S}_\mathrm{fine})-P(\boldsymbol{S}_\mathrm{coarse}), & \ell > \ell_0 \end{cases}.\]

Cancellations from the telescoping sum allow us to write

\[\lim_{\ell \to \infty} P_\ell = P_{\ell_0} + \sum_{\ell=\ell_0+1}^\infty P_\ell.\]

Examples:

>>> integrand = FinancialOption(DigitalNetB2(dimension=3,seed=7),option="EUROPEAN")
>>> y = integrand(2**10)
>>> y.shape
(1024,)
>>> print("%.4f"%y.mean())
4.2126
>>> print("%.4f"%integrand.get_exact_value())
4.2115
>>> integrand
FinancialOption (AbstractIntegrand)
    option          EUROPEAN
    call_put        CALL
    volatility      2^(-1)
    start_price     30
    strike_price    35
    interest_rate   0
    t_final         1
>>> integrand.true_measure
BrownianMotion (AbstractTrueMeasure)
    time_vec        [0.333 0.667 1.   ]
    drift           0
    mean            [0. 0. 0.]
    covariance      [[0.333 0.333 0.333]
                     [0.333 0.667 0.667]
                     [0.333 0.667 1.   ]]
    decomp_type     PCA

>>> integrand = FinancialOption(DigitalNetB2(dimension=64,seed=7),option="ASIAN")
>>> y = integrand(2**10)
>>> y.shape
(1024,)
>>> print("%.4f"%y.mean())
1.7782

With independent replications

>>> integrand = FinancialOption(DigitalNetB2(dimension=64,seed=7,replications=2**4),option="ASIAN")
>>> y = integrand(2**6)
>>> y.shape
(16, 64)
>>> muhats = y.mean(-1) 
>>> muhats.shape 
(16,)
>>> print("%.4f"%muhats.mean())
1.7765

Multi-level options

>>> seed_seq = np.random.SeedSequence(7) 
>>> d_coarsest = 8
>>> num_levels = 4
>>> ns = [2**11,2**10,2**9,2**8]
>>> integrands = [FinancialOption(DigitalNetB2(dimension=2**l*d_coarsest,seed=seed_seq.spawn(1)[0]),option="ASIAN",level=l,d_coarsest=d_coarsest) for l in range(num_levels)]
>>> ys = [integrands[l](ns[l]) for l in range(num_levels)]
>>> for l in range(num_levels):
...     print("ys[%d].shape = %s"%(l,ys[l].shape))
ys[0].shape = (2, 2048)
ys[1].shape = (2, 1024)
ys[2].shape = (2, 512)
ys[3].shape = (2, 256)
>>> ymeans = np.stack([(ys[l][1]-ys[l][0]).mean(-1) for l in range(num_levels)])
>>> with np.printoptions(formatter={"float": lambda x: "%.2e"%x}):
...     ymeans
array([1.78e+00, 2.61e-03, 5.57e-04, 1.10e-03])
>>> print("%.4f"%ymeans.sum())
1.7850

Multi-level options with independent replications

>>> seed_seq = np.random.SeedSequence(7) 
>>> d_coarsest = 8
>>> num_levels = 4
>>> ns = [2**7,2**6,2**5,2**4]
>>> integrands = [FinancialOption(DigitalNetB2(dimension=2**l*d_coarsest,seed=seed_seq.spawn(1)[0],replications=2**4),option="ASIAN",level=l,d_coarsest=d_coarsest) for l in range(num_levels)]
>>> ys = [integrands[l](ns[l]) for l in range(num_levels)]
>>> for l in range(num_levels):
...     print("ys[%d].shape = %s"%(l,ys[l].shape))
ys[0].shape = (2, 16, 128)
ys[1].shape = (2, 16, 64)
ys[2].shape = (2, 16, 32)
ys[3].shape = (2, 16, 16)
>>> muhats = np.stack([(ys[l][1]-ys[l][0]).mean(-1) for l in range(num_levels)])
>>> muhats.shape
(4, 16)
>>> muhathat = muhats.mean(-1)
>>> with np.printoptions(formatter={"float": lambda x: "%.2e"%x}):
...     muhathat
array([1.80e+00, -3.08e-03, 2.64e-03, 1.14e-03])
>>> print("%.4f"%muhathat.sum())
1.7982

References:

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

Parameters:

Name	Type	Description	Default
sampler	Union[AbstractDiscreteDistribution, AbstractTrueMeasure]	Either a discrete distribution from which to transform samples, or a true measure by which to compose a transform.	required
option	str	Option type in ['ASIAN', 'EUROPEAN', 'BARRIER', 'LOOKBACK', 'DIGITAL']	'ASIAN'
call_put	str	Either 'CALL' or 'PUT'.	'CALL'
volatility	float	\(\sigma\).	0.5
start_price	float	\(S_0\).	30
strike_price	float	\(K\).	35
interest_rate	float	\(r\).	0
t_final	float	\(\tau_d\).	1
decomp_type	str	Method for decomposition for covariance matrix. Options include 'PCA' for principal component analysis, or 'Cholesky' for cholesky decomposition.	'PCA'
level	Union[None, int]	Level for multilevel problems	None
d_coarsest	Union[None, int]	Dimension of the problem on the coarsest level.	2
asian_mean	str	Either 'ARITHMETIC' or 'GEOMETRIC'.	'ARITHMETIC'
asian_mean_quadrature_rule	str	Either 'TRAPEZOIDAL' or 'RIGHT'.	'TRAPEZOIDAL'
barrier_in_out	str	Either 'IN' or 'OUT'.	'IN'
barrier_price	float	\(B\).	38
digital_payout	float	\(\rho\).	10

Source code in qmcpy/integrand/financial_option.py

def __init__(self, 
             sampler,
             option = "ASIAN", 
             call_put = 'CALL',
             volatility = 0.5,
             start_price = 30,
             strike_price = 35,
             interest_rate = 0,
             t_final = 1,
             decomp_type = 'PCA',
             level = None,
             d_coarsest = 2,
             asian_mean = "ARITHMETIC",
             asian_mean_quadrature_rule = "TRAPEZOIDAL",
             barrier_in_out = "IN",
             barrier_price = 38,
             digital_payout = 10):
    r"""
    Args:
        sampler (Union[AbstractDiscreteDistribution,AbstractTrueMeasure]): Either  

            - a discrete distribution from which to transform samples, or
            - a true measure by which to compose a transform.
        option (str): Option type in `['ASIAN', 'EUROPEAN', 'BARRIER', 'LOOKBACK', 'DIGITAL']`
        call_put (str): Either `'CALL'` or `'PUT'`. 
        volatility (float): $\sigma$.
        start_price (float): $S_0$.
        strike_price (float): $K$.
        interest_rate (float): $r$.
        t_final (float): $\tau_d$.
        decomp_type (str): Method for decomposition for covariance matrix. Options include

            - `'PCA'` for principal component analysis, or 
            - `'Cholesky'` for cholesky decomposition.
        level (Union[None,int]): Level for multilevel problems 
        d_coarsest (Union[None,int]): Dimension of the problem on the coarsest level.
        asian_mean (str): Either `'ARITHMETIC'` or `'GEOMETRIC'`.
        asian_mean_quadrature_rule (str): Either 'TRAPEZOIDAL' or 'RIGHT'. 
        barrier_in_out (str): Either `'IN'` or `'OUT'`. 
        barrier_price (float): $B$. 
        digital_payout (float): $\rho$. 
    """
    self.parameters = ['option', 'call_put', 'volatility', 'start_price', 'strike_price', 'interest_rate', 't_final']
    self.t_final = t_final
    self.sampler = sampler
    self.decomp_type = decomp_type
    self.true_measure = BrownianMotion(self.sampler,t_final=self.t_final,decomp_type=self.decomp_type)
    self.volatility = float(volatility)
    self.start_price = float(start_price)
    self.strike_price = float(strike_price)
    self.interest_rate = float(interest_rate)
    self.discount_factor = np.exp(-self.interest_rate*self.t_final)
    self.level = level
    self.d_coarsest = d_coarsest
    if self.level is not None:
        self.multilevel = True 
        self.parameters += ['level','d_coarsest']
        assert np.isscalar(self.level) and self.level%1==0 
        assert np.isscalar(self.d_coarsest) and self.d_coarsest%1==0 and d_coarsest>0 and np.log2(d_coarsest)%1==0, "d_coarsest must be an integer power of 2"
        self.level = int(self.level)
        self.d_coarsest = int(self.d_coarsest) 
        assert self.sampler.d==self.d_coarsest*2**self.level, "the dimension of the sampler must equal d_coarsest*2^level = %d"%(d_coarsest*2**self.level)
        self.cost = self.d_coarsest*2**self.level
        dim_shape = (2,)
    else:
        self.multilevel = False
        dim_shape = ()
    self.call_put = str(call_put).upper()
    assert self.call_put in ['CALL','PUT'], "invalid call_put = %s"%self.call_put
    self.option = str(option).upper()
    self.asian_mean = str(asian_mean).upper()
    self.asian_mean_quadrature_rule = str(asian_mean_quadrature_rule).upper()
    self.barrier_in_out = str(barrier_in_out).upper()
    assert np.isscalar(barrier_price)
    self.barrier_price = float(barrier_price)
    assert np.isscalar(digital_payout) and digital_payout>0
    self.digital_payout = float(digital_payout)
    if self.option=="EUROPEAN":
        self.payoff = self.payoff_european_call if self.call_put=='CALL' else self.payoff_european_put
    elif self.option=="ASIAN":
        self.parameters += ['asian_mean']
        assert self.asian_mean in ['ARITHMETIC','GEOMETRIC'], "invalid asian_mean = %s"%self.asian_mean
        assert self.asian_mean_quadrature_rule in ['TRAPEZOIDAL','RIGHT'], "invalid asian_mean_quadrature_rule = %s"%self.asian_mean_quadrature_rule
        if self.asian_mean=="ARITHMETIC":
            if self.asian_mean_quadrature_rule=="TRAPEZOIDAL":
                self.payoff = self.payoff_asian_arithmetic_trap_call if self.call_put=='CALL' else self.payoff_asian_arithmetic_trap_put
            elif self.asian_mean_quadrature_rule=="RIGHT":
                self.payoff = self.payoff_asian_arithmetic_right_call if self.call_put=='CALL' else self.payoff_asian_arithmetic_right_put
        elif self.asian_mean=="GEOMETRIC":
            if self.asian_mean_quadrature_rule=="TRAPEZOIDAL":
                self.payoff = self.payoff_asian_geometric_trap_call if self.call_put=='CALL' else self.payoff_asian_geometric_trap_put
            elif self.asian_mean_quadrature_rule=="RIGHT":
                self.payoff = self.payoff_asian_geometric_right_call if self.call_put=='CALL' else self.payoff_asian_geometric_right_put
    elif self.option=="BARRIER":
        self.parameters += ['barrier_in_out','barrier_up_down','barrier_price']
        self.barrier_up_down = "UP" if self.start_price<self.barrier_price else "DOWN"
        if self.barrier_in_out=="IN":
            if self.barrier_up_down=="UP":
                self.payoff = self.payoff_barrier_in_up_call if self.call_put=='CALL' else self.payoff_barrier_in_up_put
            else:
                self.payoff = self.payoff_barrier_in_down_call if self.call_put=='CALL' else self.payoff_barrier_in_down_put
        elif self.barrier_in_out=="OUT":
            if self.barrier_up_down=="UP":
                self.payoff = self.payoff_barrier_out_up_call if self.call_put=='CALL' else self.payoff_barrier_out_up_put
            else:
                self.payoff = self.payoff_barrier_out_down_call if self.call_put=='CALL' else self.payoff_barrier_out_down_put
        else:
            raise ParameterError("invalid barrier_in_out = %s"%self.barrier_in_out)
    elif self.option=="LOOKBACK":
        self.payoff = self.payoff_lookback_call if self.call_put=='CALL' else self.payoff_lookback_put
    elif self.option=="DIGITAL":
        self.payoff = self.payoff_digital_call if self.call_put=='CALL' else self.payoff_digital_put
    else:
        raise ParameterError("invalid option type %s"%self.option)
    super(FinancialOption,self).__init__(dimension_indv=dim_shape,dimension_comb=dim_shape,parallel=False)  

get_exact_value

get_exact_value()

Compute the exact analytic fair price of the option in finite dimensions. Supports

• option='EUROPEAN'
• option='ASIAN' with asian_mean='GEOMETRIC' and asian_mean_quadrature_rule='RIGHT'

Returns:

Name	Type	Description
mean	float	Exact value of the integral.

Source code in qmcpy/integrand/financial_option.py

def get_exact_value(self):
    """
    Compute the exact analytic fair price of the option in finite dimensions. Supports 

    - `option='EUROPEAN'`
    - `option='ASIAN'` with `asian_mean='GEOMETRIC'` and `asian_mean_quadrature_rule='RIGHT'`

    Returns: 
        mean (float): Exact value of the integral. 
    """
    if self.option=="EUROPEAN":
        denom = self.volatility*np.sqrt(self.t_final)
        decay = self.strike_price*self.discount_factor
        if self.call_put == 'CALL':
            term1 = np.log(self.start_price/self.strike_price)+(self.interest_rate+self.volatility**2/2)*self.t_final
            term2 = np.log(self.start_price/self.strike_price)+(self.interest_rate-self.volatility**2/2)*self.t_final
            fp = self.start_price*norm.cdf(term1/denom)-decay*norm.cdf(term2/denom)
        elif self.call_put == 'PUT':
            term1 = np.log(self.strike_price/self.start_price)-(self.interest_rate-self.volatility**2/2)*self.t_final
            term2 = np.log(self.strike_price/self.start_price)-(self.interest_rate+self.volatility**2/2)*self.t_final
            fp = decay*norm.cdf(term1/denom)-self.start_price*norm.cdf(term2/denom)
    elif self.option=="ASIAN":
        assert self.asian_mean=='GEOMETRIC' and self.asian_mean_quadrature_rule=='RIGHT', "exact value for Asian options only implemented for self.asian_mean=='GEOMETRIC' and self.asian_mean_quadrature_rule=='RIGHT'"
        Tbar = (1+1/self.d)*self.t_final/2
        sigmabar = self.volatility*np.sqrt((2+1/self.d)/3)
        rbar = self.interest_rate+(sigmabar**2-self.volatility**2)/2
        gmeancall,gmeanput = _eurogbmprice(self.start_price,rbar,Tbar,sigmabar,self.strike_price)
        if self.call_put=='CALL':
            fp = gmeancall * np.exp(rbar*Tbar-self.interest_rate*self.t_final)
        elif self.call_put=='PUT':
            fp = gmeanput * np.exp(rbar*Tbar-self.interest_rate*self.t_final)
        return fp
    else:
        raise ParameterError("exact value not supported for option = %s"%self.option)
    return fp

get_exact_value_inf_dim

get_exact_value_inf_dim()

Get the exact analytic fair price of the option in infinite dimensions. Supports

• option='ASIAN' with asian_mean='GEOMETRIC'

Returns:

Name	Type	Description
mean	float	Exact value of the integral.

Source code in qmcpy/integrand/financial_option.py

def get_exact_value_inf_dim(self):
    r"""
    Get the exact analytic fair price of the option in infinite dimensions. Supports 

    - `option='ASIAN'` with `asian_mean='GEOMETRIC'`

    Returns: 
        mean (float): Exact value of the integral. 
    """
    if self.option=='ASIAN':
        assert self.asian_mean=='GEOMETRIC' , "get_exact_value_inf_dim for the Asian option only available for self.asian_mean=='GEOMETRIC'"
        sigma_g = self.volatility/np.sqrt(3) 
        b = 1/2*(self.interest_rate-1/2*sigma_g**2)
        d1 = (np.log(self.start_price/self.strike_price)+(b+1/2*sigma_g**2)*self.t_final)/(sigma_g*np.sqrt(self.t_final))
        d2 = d1-sigma_g*np.sqrt(self.t_final)
        f1 = self.start_price*np.exp((b-self.interest_rate)*self.t_final)
        f2 = self.strike_price*np.exp(-self.interest_rate*self.t_final)
        if self.call_put=="CALL":
            val = f1*norm.cdf(d1)-f2*norm.cdf(d2)
        elif self.call_put=="PUT":
            val = f2*norm.cdf(-d2)-f1*norm.cdf(-d1)
    else:
        raise Exception("get_exact_value_inf_dim not implemented for option = %s"%self.option)
    return val

CustomFun

Bases: AbstractIntegrand

User supplied integrand \(g\). In the following example we implement

Examples:

First we will implement

\[g(\boldsymbol{t}) = t_1^2t_2, \qquad \boldsymbol{T}=(T_1,T_2) \sim \mathcal{N}((1,2)^T,\mathsf{I}).\]

>>> integrand = CustomFun(
...     true_measure = Gaussian(DigitalNetB2(2,seed=7),mean=[1,2]),
...     g = lambda t: t[...,0]**2*t[...,1])
>>> y = integrand(2**10)
>>> print("%.4f"%y.mean())
3.9991

With independent replications

>>> integrand = CustomFun(
...     true_measure = Gaussian(DigitalNetB2(2,seed=7,replications=2**4),mean=[1,2]),
...     g = lambda t: t[...,0]**2*t[...,1])
>>> y = integrand(2**6)
>>> y.shape
(16, 64)
>>> muhats = y.mean(1) 
>>> muhats.shape 
(16,)
>>> print("%.4f"%muhats.mean())
3.9330

Next we will implement the multi-output function

\[g(\boldsymbol{t}) = \begin{pmatrix} \sin(t_1)\cos(t_2) \\ \cos(t_1)\sin(t_2) \\ \sin(t_1)+\cos(t_2) \\ \cos(t_1)+\sin(t_2) \end{pmatrix} \qquad \boldsymbol{T}=(T_1,T_2) \sim \mathcal{U}[0,2\pi]^2.\]

>>> def g(t):
...     t1,t2 = t[...,0],t[...,1]
...     sint1,cost1,sint2,cost2 = np.sin(t1),np.cos(t1),np.sin(t2),np.cos(t2)
...     y1 = sint1*cost2
...     y2 = cost1*sint2
...     y3 = sint1+cost2
...     y4 = cost1+sint2
...     y = np.stack([y1,y2,y3,y4])
...     return y
>>> integrand = CustomFun(
...     true_measure = Uniform(DigitalNetB2(2,seed=7),lower_bound=0,upper_bound=2*np.pi),
...     g = g, 
...     dimension_indv = (4,))
>>> x = integrand.discrete_distrib(2**10)
>>> y = integrand.f(x)
>>> y.shape
(4, 1024)
>>> with np.printoptions(formatter={"float": lambda x: "%.2e"%x}):
...     y.mean(-1)
array([8.18e-04, 1.92e-06, -2.26e-10, 5.05e-07])

Stopping criterion which supporting vectorized outputs may pass in Boolean compute_flags with dimension_indv shape indicating which output need to evaluated, i.e. where compute_flags is False we do not need to evaluate the integrand. We have not used this in inexpensive example above.

With independent replications

>>> integrand = CustomFun(
...     true_measure = Uniform(DigitalNetB2(2,seed=7,replications=2**4),lower_bound=0,upper_bound=2*np.pi),
...     g = g, 
...     dimension_indv = (4,))
>>> x = integrand.discrete_distrib(2**6)
>>> x.shape
(16, 64, 2)
>>> y = integrand.f(x)
>>> y.shape
(4, 16, 64)
>>> muhats = y.mean(-1) 
>>> muhats.shape 
(4, 16)
>>> with np.printoptions(formatter={"float": lambda x: "%.2e"%x}):
...     muhats.mean(-1)
array([3.83e-03, -6.78e-03, -1.56e-03, -5.65e-04])

Parameters:

Name	Type	Description	Default
true_measure	AbstractTrueMeasure	The true measure.	required
g	callable	A function handle.	required
dimension_indv	tuple	Shape of individual solution outputs from g.	()
parallel	int	Parallelization flag. When parallel = 0 or parallel = 1 then function evaluation is done in serial fashion. parallel > 1 specifies the number of processes used by multiprocessing.Pool or multiprocessing.pool.ThreadPool. Setting parallel=True is equivalent to parallel = os.cpu_count().	False

Note

For parallel > 1 do not set g to be anonymous function (i.e. a lambda function)

Source code in qmcpy/integrand/custom_fun.py

def __init__(self, true_measure, g, dimension_indv=(), parallel=False):
    """
    Args:
        true_measure (AbstractTrueMeasure): The true measure. 
        g (callable): A function handle. 
        dimension_indv (tuple): Shape of individual solution outputs from `g`.
        parallel (int): Parallelization flag. 

            - When `parallel = 0` or `parallel = 1` then function evaluation is done in serial fashion.
            - `parallel > 1` specifies the number of processes used by `multiprocessing.Pool` or `multiprocessing.pool.ThreadPool`.

            Setting `parallel=True` is equivalent to `parallel = os.cpu_count()`.

    Note:
        For `parallel > 1` do *not* set `g` to be anonymous function (i.e. a `lambda` function)
    """
    self.parameters = []
    self.true_measure = true_measure
    self.sampler = self.true_measure
    self.__g = g            
    super(CustomFun,self).__init__(dimension_indv=dimension_indv,dimension_comb=dimension_indv,parallel=parallel)

UMBridgeWrapper

Bases: AbstractIntegrand

Wrapper around a UM-Bridge model. See also the UM-Bridge documentation for the QMCPy client. Requires Docker is installed.

Examples:

>>> _ = os.system('docker run --name muqbppytest -dit -p 4243:4243 linusseelinger/benchmark-muq-beam-propagation:latest > /dev/null')
>>> import umbridge
>>> dnb2 = DigitalNetB2(dimension=3,seed=7)
>>> true_measure = Uniform(dnb2,lower_bound=1,upper_bound=1.05)
>>> um_bridge_model = umbridge.HTTPModel('http://localhost:4243','forward')
>>> um_bridge_config = {"d": dnb2.d}
>>> integrand = UMBridgeWrapper(true_measure,um_bridge_model,um_bridge_config,parallel=False)
>>> y = integrand(2**10)
>>> with np.printoptions(formatter={"float":lambda x: "%.1e"%x}):
...     y.mean(-1)
array([0.0e+00, 3.9e+00, 1.5e+01, 3.2e+01, 5.5e+01, 8.3e+01, 1.2e+02,
       1.5e+02, 2.0e+02, 2.4e+02, 2.9e+02, 3.4e+02, 3.9e+02, 4.3e+02,
       4.7e+02, 5.0e+02, 5.3e+02, 5.6e+02, 5.9e+02, 6.2e+02, 6.4e+02,
       6.6e+02, 6.9e+02, 7.2e+02, 7.6e+02, 7.9e+02, 8.3e+02, 8.6e+02,
       9.0e+02, 9.4e+02, 9.7e+02])
>>> _ = os.system('docker rm -f muqbppytest > /dev/null')

Custom model with independent replications

>>> class TestModel(umbridge.Model):
...     def __init__(self):
...         super().__init__("forward")
...     def get_input_sizes(self, config):
...         return [1,2,3]
...     def get_output_sizes(self, config):
...         return [3,2,1]
...     def __call__(self, parameters, config):
...         out0 = [parameters[2][0],sum(parameters[2][:2]),sum(parameters[2])]
...         out1 = [parameters[1][0],sum(parameters[1])]
...         out2 = [parameters[0]]
...         return [out0,out1,out2]
...     def supports_evaluate(self):
...         return True
>>> um_bridge_model = TestModel()
>>> um_bridge_config = {}
>>> d = sum(um_bridge_model.get_input_sizes(config=um_bridge_config))
>>> true_measure = Uniform(DigitalNetB2(dimension=d,seed=7,replications=15),lower_bound=-1,upper_bound=1)
>>> integrand = UMBridgeWrapper(true_measure,um_bridge_model)
>>> y = integrand(2**6)
>>> y.shape 
(6, 15, 64)
>>> muhats = y.mean(-1)
>>> muhats.shape 
(6, 15)
>>> muhats_aggregate = muhats.mean(-1)
>>> muhats_aggregate.shape 
(6,)
>>> muhats_agg_list_of_lists = integrand.to_umbridge_out_sizes(muhats_aggregate)
>>> [["%.2e"%ii for ii in i] for i in muhats_agg_list_of_lists]
[['-1.59e-08', '1.49e-04', '1.49e-04'], ['8.20e-06', '-1.38e-04'], ['-8.14e-06']]

Parameters:

Name	Type	Description	Default
true_measure	AbstractTrueMeasure	The true measure.	required
model	HTTPModel	A UM-Bridge model.	required
config	dict	Configuration keyword argument to umbridge.HTTPModel(url,name).__call__.	{}
parallel	int	Parallelization flag. When parallel = 0 or parallel = 1 then function evaluation is done in serial fashion. parallel > 1 specifies the number of processes used by multiprocessing.Pool or multiprocessing.pool.ThreadPool. Setting parallel=True is equivalent to parallel = os.cpu_count().	False

Source code in qmcpy/integrand/umbridge_wrapper.py

def __init__(self, true_measure, model, config={}, parallel=False):
    """
    Args:
        true_measure (AbstractTrueMeasure): The true measure.  
        model (umbridge.HTTPModel): A `UM-Bridge` model. 
        config (dict): Configuration keyword argument to `umbridge.HTTPModel(url,name).__call__`.
        parallel (int): Parallelization flag. 

            - When `parallel = 0` or `parallel = 1` then function evaluation is done in serial fashion.
            - `parallel > 1` specifies the number of processes used by `multiprocessing.Pool` or `multiprocessing.pool.ThreadPool`.

            Setting `parallel=True` is equivalent to `parallel = os.cpu_count()`.
    """
    import umbridge
    self.parameters = []
    self.true_measure = true_measure
    self.sampler = self.true_measure 
    self.model = model
    if not self.model.supports_evaluate(): raise ParameterError("UMBridgeWrapper requires model supports evaluation.")
    self.config = config
    self.parallel = parallel
    self.d_in_umbridge = np.append(0,np.cumsum(self.model.get_input_sizes(self.config)))
    self.n_d_in_umbridge = len(self.d_in_umbridge)-1
    if self.true_measure.d!=self.d_in_umbridge[-1]:
        raise ParameterError("sampler dimension (%d) must equal the sum of UMBridgeWrapper input sizes (%d)."%(self.true_measure.d,self.d_in_umbridge[-1]))
    self.d_out_umbridge = np.append(0,np.cumsum(self.model.get_output_sizes(self.config)))
    self.n_d_out_umbridge = len(self.d_out_umbridge)-1
    self.total_out_elements = int(self.d_out_umbridge[-1])
    super(UMBridgeWrapper,self).__init__(
        dimension_indv = () if self.total_out_elements==1 else (self.total_out_elements,),
        dimension_comb = () if self.total_out_elements==1 else (self.total_out_elements,),
        parallel = self.parallel,
        threadpool = True)

to_umbridge_out_sizes

to_umbridge_out_sizes(x)

Convert a data attribute to UM-Bridge output sized list of lists.

Parameters:

Name	Type	Description	Default
x	ndarray	Array of length sum(model.get_output_sizes(self.config)) where model is a umbridge.HTTPModel.	required

Returns:

Name	Type	Description
x_list_list	list	List of lists with sub-list lengths specified by model.get_output_sizes(self.config).

Source code in qmcpy/integrand/umbridge_wrapper.py

def to_umbridge_out_sizes(self, x):
    """
    Convert a data attribute to `UM-Bridge` output sized list of lists. 

    Args:
        x (np.ndarray): Array of length `sum(model.get_output_sizes(self.config))` where `model` is a `umbridge.HTTPModel`. 

    Returns:
        x_list_list (list): List of lists with sub-list lengths specified by `model.get_output_sizes(self.config)`.
    """
    return [x[...,self.d_out_umbridge[j]:self.d_out_umbridge[j+1]].tolist() for j in range(self.n_d_out_umbridge)]

SensitivityIndices

Bases: AbstractIntegrand

Sensitivity indices i.e. normalized Sobol' Indices.

Examples:

Singleton indices

>>> function = Keister(DigitalNetB2(dimension=4,seed=7))
>>> integrand = SensitivityIndices(function,indices='singletons')
>>> integrand.indices
array([[ True, False, False, False],
       [False,  True, False, False],
       [False, False,  True, False],
       [False, False, False,  True]])
>>> y = integrand(2**10)
>>> y.shape
(2, 3, 4, 1024)
>>> ymean = y.mean(-1)
>>> ymean.shape
(2, 3, 4)
>>> sigma_hat = ymean[:,2,:]-ymean[:,1,:]**2
>>> sigma_hat.shape 
(2, 4)
>>> sigma_hat
array([[17.81421941, 17.81421941, 17.81421941, 17.81421941],
       [17.81421941, 17.81421941, 17.81421941, 17.81421941]])
>>> closed_total_approx = ymean[:,0]/sigma_hat
>>> closed_total_approx.shape 
(2, 4)
>>> closed_total_approx
array([[0.23756849, 0.2423556 , 0.24844127, 0.24771739],
       [0.25090825, 0.25344588, 0.24531633, 0.26244379]])

Check what all indices look like for \(d=3\)

>>> integrand = SensitivityIndices(Keister(DigitalNetB2(dimension=3,seed=7)),indices='all')
>>> integrand.indices
array([[ True, False, False],
       [False,  True, False],
       [False, False,  True],
       [ True,  True, False],
       [ True, False,  True],
       [False,  True,  True]])

Vectorized function for all singletons and pairs of dimensions

>>> function = BoxIntegral(DigitalNetB2(dimension=4,seed=7,replications=2**4),s=np.arange(1,31).reshape((5,6)))
>>> indices = np.zeros((function.d,function.d,function.d),dtype=bool) 
>>> r = np.arange(function.d) 
>>> indices[r,:,r] = True 
>>> indices[:,r,r] = True 
>>> integrand = SensitivityIndices(function,indices=indices)
>>> integrand.indices.shape
(4, 4, 4)
>>> integrand.indices
array([[[ True, False, False, False],
        [ True,  True, False, False],
        [ True, False,  True, False],
        [ True, False, False,  True]],

       [[ True,  True, False, False],
        [False,  True, False, False],
        [False,  True,  True, False],
        [False,  True, False,  True]],

       [[ True, False,  True, False],
        [False,  True,  True, False],
        [False, False,  True, False],
        [False, False,  True,  True]],

       [[ True, False, False,  True],
        [False,  True, False,  True],
        [False, False,  True,  True],
        [False, False, False,  True]]])
>>> y = integrand(2**10)
>>> y.shape 
(2, 3, 4, 4, 5, 6, 16, 1024)
>>> muhats = y.mean(-1) 
>>> muhats.shape 
(2, 3, 4, 4, 5, 6, 16)
>>> muhathat = muhats.mean(-1) 
>>> muhathat.shape 
(2, 3, 4, 4, 5, 6)
>>> sigma_hat = muhathat[:,2,:]-muhathat[:,1,:]**2
>>> sigma_hat.shape
(2, 4, 4, 5, 6)
>>> closed_total_approx = muhathat[:,0]/sigma_hat
>>> closed_total_approx.shape
(2, 4, 4, 5, 6)

References:

1. Aleksei G. Sorokin and Jagadeeswaran Rathinavel.
   On Bounding and Approximating Functions of Multiple Expectations Using Quasi-Monte Carlo.
   International Conference on Monte Carlo and Quasi-Monte Carlo Methods in Scientific Computing.
   Cham: Springer International Publishing, 2022.
   https://link.springer.com/chapter/10.1007/978-3-031-59762-6_29.

2. Art B. Owen.
   Monte Carlo theory, methods and examples.
   Appendix A. Equations (A.16) and (A.18). 2013. https://artowen.su.domains/mc/A-anova.pdf.

Parameters:

Name	Type	Description	Default
integrand	AbstractIntegrand	Integrand to find sensitivity indices of.	required
indices	ndarray	Bool array with shape \((\dots,d)\) where each length \(d\) vector item indicates which dimensions are active in the subset. The default indices='singletons' sets indices=np.eye(d,dtype=bool). Setting incides='all' sets indices = np.array([[bool(int(b)) for b in np.binary_repr(i,width=d)] for i in range(1,2**d-1)],dtype=bool)	'singletons'

Source code in qmcpy/integrand/sensitivity_indices.py

def __init__(self, integrand, indices='singletons'):
    r"""
    Args:
        integrand (AbstractIntegrand): Integrand to find sensitivity indices of.
        indices (np.ndarray): Bool array with shape $(\dots,d)$ where each length $d$ vector item indicates which dimensions are active in the subset.

            - The default `indices='singletons'` sets `indices=np.eye(d,dtype=bool)`.
            - Setting `incides='all'` sets `indices = np.array([[bool(int(b)) for b in np.binary_repr(i,width=d)] for i in range(1,2**d-1)],dtype=bool)`
    """
    self.parameters = ['indices']
    self.integrand = integrand
    self.dtilde = self.integrand.d
    assert self.dtilde>1, "SensitivityIndices does not make sense for d=1"
    self.indices = indices
    if isinstance(self.indices,str) and self.indices=='singletons':
        self.indices = np.eye(self.dtilde,dtype=bool)
    elif isinstance(self.indices,str) and self.indices=='all':
        self.indices = np.zeros((0,self.dtilde),dtype=bool)
        for r in range(1,self.dtilde):
            idxs_r = np.zeros((int(scipy.special.comb(self.dtilde,r)),self.dtilde),dtype=bool)
            for i,comb in enumerate(combinations(range(self.dtilde),r)):
                idxs_r[i,comb] = True 
            self.indices = np.vstack([self.indices,idxs_r])
    self.indices = np.atleast_1d(self.indices)
    assert self.indices.dtype==bool and self.indices.ndim>=1 and self.indices.shape[-1]==self.dtilde 
    assert not (self.indices==self.indices[...,0,None]).all(-1).any(), "indices cannot include the emptyset or the set of all dimensions"
    self.not_indices = ~self.indices
    # sensitivity_index
    self.true_measure = self.integrand.true_measure
    self.discrete_distrib = self.integrand.discrete_distrib.spawn(s=1,dimensions=[2*self.dtilde])[0]
    self.sampler = self.integrand.sampler
    self.i_slice = (slice(None),)*len(self.integrand.d_indv)
    super(SensitivityIndices,self).__init__(
        dimension_indv = (2,3)+self.indices.shape[:-1]+self.integrand.d_indv,
        dimension_comb = (2,)+self.indices.shape[:-1]+self.integrand.d_indv,
        parallel = False)
    self.d = 2*self.dtilde

Keister

Bases: AbstractIntegrand

Keister function from [1].

\[f(\boldsymbol{t}) = \pi^{d/2} \cos(\lVert \boldsymbol{t} \rVert_2) \qquad \boldsymbol{T} \sim \mathcal{N}(\boldsymbol{0},\mathsf{I}/2).\]

Examples:

>>> integrand = Keister(DigitalNetB2(2,seed=7))
>>> y = integrand(2**10)
>>> print("%.4f"%y.mean())
1.8080
>>> integrand.true_measure
Gaussian (AbstractTrueMeasure)
    mean            0
    covariance      2^(-1)
    decomp_type     PCA

With independent replications

>>> integrand = Keister(DigitalNetB2(2,seed=7,replications=2**4))
>>> y = integrand(2**6)
>>> y.shape
(16, 64)
>>> muhats = y.mean(-1) 
>>> muhats.shape 
(16,)
>>> print("%.4f"%muhats.mean())
1.8024

References:

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

Parameters:

Name	Type	Description	Default
sampler	Union[AbstractDiscreteDistribution, AbstractTrueMeasure]	Either a discrete distribution from which to transform samples, or a true measure by which to compose a transform.	required

Source code in qmcpy/integrand/keister.py

def __init__(self, sampler):
    r"""
    Args:
        sampler (Union[AbstractDiscreteDistribution,AbstractTrueMeasure]): Either  

            - a discrete distribution from which to transform samples, or
            - a true measure by which to compose a transform.
    """
    self.sampler = sampler
    self.true_measure = Gaussian(self.sampler,mean=0,covariance=1/2)
    super(Keister,self).__init__(dimension_indv=(),dimension_comb=(),parallel=False)

get_exact_value classmethod

get_exact_value(d)

Compute the exact analytic value of the Keister integral with dimension \(d\).

Parameters:

Name	Type	Description	Default
d	int	Dimension.	required

Returns:

Name	Type	Description
mean	float	Exact value of the integral.

Source code in qmcpy/integrand/keister.py

@classmethod
def get_exact_value(self, d):
    """
    Compute the exact analytic value of the Keister integral with dimension $d$. 

    Args:
        d (int): Dimension. 

    Returns: 
        mean (float): Exact value of the integral. 
    """
    cosinteg = np.zeros(shape=(d))
    cosinteg[0] = np.sqrt(np.pi) / (2 * np.exp(1 / 4))
    sininteg = np.zeros(shape=(d))
    sininteg[0] = 4.244363835020225e-01
    cosinteg[1] = (1 - sininteg[0]) / 2
    sininteg[1] = cosinteg[0] / 2
    for j in range(2, d):
        cosinteg[j] = ((j-1)*cosinteg[j-2]-sininteg[j-1])/2
        sininteg[j] = ((j-1)*sininteg[j-2]+cosinteg[j-1])/2
    I = (2*(np.pi**(d/2))/gamma(d/2))*cosinteg[d-1]
    return I

Genz

Bases: AbstractIntegrand

Genz function following the DAKOTA implementation.

\[g_\mathrm{oscillatory}(\boldsymbol{t}) = \cos\left(-\sum_{j=1}^d c_j t_j\right)\]

or

\[g_\mathrm{corner-peak}(\boldsymbol{t}) = \left(1+\sum_{j=1}^d c_j t_j\right)^{-(d+1)}\]

where

\[\boldsymbol{T} \sim \mathcal{U}[0,1]^d\]

and the coefficients \(\boldsymbol{c}\) are have three kinds

\[c_k^{(1)} = \frac{k-1/2}{d}, \qquad c_k^{(2)} = \frac{1}{k^2}, \qquad c_k^{(3)} = \exp\left(\frac{k \log(10^{-8})}{d}\right), \qquad k=1,\dots,d.\]

Examples:

>>> for kind_func in ['OSCILLATORY','CORNER PEAK']:
...     for kind_coeff in [1,2,3]:
...         integrand = Genz(DigitalNetB2(2,seed=7),kind_func=kind_func,kind_coeff=kind_coeff)
...         y = integrand(2**14)
...         mu_hat = y.mean()
...         print('%-15s %-3d %.3f'%(kind_func,kind_coeff,mu_hat))
OSCILLATORY     1   -0.351
OSCILLATORY     2   -0.329
OSCILLATORY     3   -0.217
CORNER PEAK     1   0.713
CORNER PEAK     2   0.714
CORNER PEAK     3   0.720

With independent replications

>>> integrand = Genz(DigitalNetB2(2,seed=7,replications=2**4),kind_func="CORNER PEAK",kind_coeff=3)
>>> y = integrand(2**6)
>>> y.shape
(16, 64)
>>> muhats = y.mean(-1) 
>>> muhats.shape 
(16,)
>>> print("%.4f"%muhats.mean())
0.7200

Parameters:

Name	Type	Description	Default
sampler	Union[AbstractDiscreteDistribution, AbstractTrueMeasure]	Either a discrete distribution from which to transform samples, or a true measure by which to compose a transform.	required
kind_func	str	Either 'OSCILLATORY' or 'CORNER PEAK'	'OSCILLATORY'
kind_coeff	int	1, 2, or 3 for choice of coefficients	1

Source code in qmcpy/integrand/genz.py

def __init__(self, sampler, kind_func='OSCILLATORY', kind_coeff=1):
    """
    Args:
        sampler (Union[AbstractDiscreteDistribution,AbstractTrueMeasure]): Either  

            - a discrete distribution from which to transform samples, or
            - a true measure by which to compose a transform.
        kind_func (str): Either `'OSCILLATORY'` or `'CORNER PEAK'`
        kind_coeff (int): 1, 2, or 3 for choice of coefficients 
    """
    self.kind_func = str(kind_func).upper().strip().replace("_"," ").replace("-"," ")
    self.kind_coeff = kind_coeff
    if (self.kind_func not in ['OSCILLATORY','CORNER PEAK']) or (self.kind_coeff not in [1,2,3]):
        raise ParameterError('''
            Genz expects 
                kind_func in ['OSCILLATORY','CORNER PEAK'] and 
                kind_coeffs in [1,2,3]''')
    self.sampler = sampler
    self.true_measure = Uniform(self.sampler)
    self.d = self.true_measure.d
    if self.kind_coeff==1: self.c = (np.arange(1,self.d+1)-.5)/self.d
    elif self.kind_coeff==2: self.c = 1/np.arange(1,self.d+1)**2
    elif self.kind_coeff==3: self.c = np.exp(np.arange(1,self.d+1)*np.log(10**(-8))/self.d)
    if self.kind_func=='OSCILLATORY':
        self.g = self.g_oscillatory
        self.c = 4.5*self.c/self.c.sum()
    elif self.kind_func=='CORNER PEAK':
        self.g = self.g_corner_peak
        self.c = 0.25*self.c/self.c.sum()
    self.c = self.c[None,:]
    self.parameters = ['kind_func','kind_coeff']
    super(Genz,self).__init__(dimension_indv=(),dimension_comb=(),parallel=False)

BoxIntegral

Bases: AbstractIntegrand

Box integral from [1], see also

\[B_s(\boldsymbol{t}) = \left(\sum_{j=1}^d t_j^2 \right)^{s/2}, \qquad \boldsymbol{T} \sim \mathcal{U}[0,1]^d.\]

Examples:

Scalar s

>>> integrand = BoxIntegral(DigitalNetB2(2,seed=7),s=7)
>>> y = integrand(2**10)
>>> y.shape
(1024,)
>>> print("%.4f"%y.mean(0))
0.7519

With independent replications

>>> integrand = BoxIntegral(DigitalNetB2(2,seed=7,replications=2**4),s=7)
>>> y = integrand(2**10)
>>> y.shape
(16, 1024)
>>> muhats = y.mean(1) 
>>> muhats.shape 
(16,)
>>> print("%.4f"%muhats.mean(0))
0.7518

Array s

>>> integrand = BoxIntegral(DigitalNetB2(5,seed=7),s=np.arange(6).reshape((2,3)))
>>> y = integrand(2**10)
>>> y.shape
(2, 3, 1024)
>>> y.mean(-1)
array([[1.        , 1.26234461, 1.66666661],
       [2.28201516, 3.22195096, 4.67188113]])

With independent replications

>>> integrand = BoxIntegral(DigitalNetB2(2,seed=7,replications=2**4),s=np.arange(6).reshape((2,3)))
>>> y = integrand(2**10)
>>> y.shape
(2, 3, 16, 1024)
>>> muhats = y.mean(-1) 
>>> muhats.shape 
(2, 3, 16)
>>> muhats.mean(-1)
array([[1.        , 0.76519118, 0.66666666],
       [0.62718785, 0.62224086, 0.64273341]])

References:

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

Parameters:

Name	Type	Description	Default
sampler	Union[AbstractDiscreteDistribution, AbstractTrueMeasure]	Either a discrete distribution from which to transform samples, or a true measure by which to compose a transform.	required
s	Union[float, ndarray]	s parameter or parameters. The output shape of g is the shape of s.	1

Source code in qmcpy/integrand/box_integral.py

def __init__(self, sampler, s=1):
    r"""
    Args:
        sampler (Union[AbstractDiscreteDistribution,AbstractTrueMeasure]): Either  

            - a discrete distribution from which to transform samples, or
            - a true measure by which to compose a transform.
        s (Union[float,np.ndarray]): `s` parameter or parameters. The output shape of `g` is the shape of `s`. 
    """
    self.parameters = ['s']
    self.s = np.array(s)
    assert self.s.size>0
    self.sampler = sampler
    self.true_measure = Uniform(self.sampler)
    self.s_over_2 = self.s/2
    super(BoxIntegral,self).__init__(dimension_indv=self.s.shape,dimension_comb=self.s.shape,parallel=False)

Ishigami

Bases: AbstractIntegrand

Ishigami function in \(d=3\) dimensions from [1] and https://www.sfu.ca/~ssurjano/ishigami.html.

\[g(\boldsymbol{t}) = (1+bt_2^4)\sin(t_0)+a\sin^2(t_1), \qquad \boldsymbol{T} = (T_0,T_1,T_2) \sim \mathcal{U}(-\pi,\pi)^3.\]

Examples:

>>> integrand = Ishigami(DigitalNetB2(3,seed=7))
>>> y = integrand(2**10)
>>> y.shape 
(1024,)
>>> print("%.4f"%y.mean())
3.5000
>>> integrand.true_measure
Uniform (AbstractTrueMeasure)
    lower_bound     -3.142
    upper_bound     3.142

With independent replications

>>> integrand = Ishigami(DigitalNetB2(3,seed=7,replications=2**4))
>>> y = integrand(2**6)
>>> y.shape
(16, 64)
>>> muhats = y.mean(-1) 
>>> muhats.shape 
(16,)
>>> print("%.4f"%muhats.mean())
3.4646

References:

1. Ishigami, T., & Homma, T.
   An importance quantification technique in uncertainty analysis for computer models.
   In Uncertainty Modeling and Analysis, 1990.
   Proceedings, First International Symposium on (pp. 398-403). IEEE.

Parameters:

Name	Type	Description	Default
sampler	Union[AbstractDiscreteDistribution, AbstractTrueMeasure]	Either a discrete distribution from which to transform samples, or a true measure by which to compose a transform.	required
a	float	First parameter \(a\).	7
b	float	Second parameter \(b\).	0.1

Source code in qmcpy/integrand/ishigami.py

def __init__(self,sampler, a=7, b=.1):
    r"""
    Args:
        sampler (Union[AbstractDiscreteDistribution,AbstractTrueMeasure]): Either  

            - a discrete distribution from which to transform samples, or
            - a true measure by which to compose a transform.
        a (float): First parameter $a$.
        b (float): Second parameter $b$.
    """
    self.sampler = sampler
    if self.sampler.d != 3:
        raise ParameterError("Ishigami integrand requires 3 dimensional sampler")
    self.a = a
    self.b = b
    self.true_measure = Uniform(self.sampler, lower_bound=-np.pi, upper_bound=np.pi)
    super(Ishigami,self).__init__(dimension_indv=(),dimension_comb=(),parallel=False)

Multimodal2d

Bases: AbstractIntegrand

Multimodal function in \(d=2\) dimensions.

\[g(\boldsymbol{t}) = (t_0^2+4)(t_1-1)/20-\sin(5t_0/2)-2 \qquad \boldsymbol{T} = (T_0,T_1) \sim \mathcal{U}([-4,7] \times [-3,8]).\]

Examples:

>>> integrand = Multimodal2d(DigitalNetB2(2,seed=7))
>>> y = integrand(2**10)
>>> print("%.4f"%y.mean())
-0.7365
>>> integrand.true_measure
Uniform (AbstractTrueMeasure)
    lower_bound     [-4 -3]
    upper_bound     [7 8]

With independent replications

>>> integrand = Multimodal2d(DigitalNetB2(2,seed=7,replications=2**4))
>>> y = integrand(2**6)
>>> y.shape
(16, 64)
>>> muhats = y.mean(-1) 
>>> muhats.shape 
(16,)
>>> print("%.4f"%muhats.mean())
-0.7366

Parameters:

Name	Type	Description	Default
sampler	Union[AbstractDiscreteDistribution, AbstractTrueMeasure]	Either a discrete distribution from which to transform samples, or a true measure by which to compose a transform.	required

Source code in qmcpy/integrand/multimodal2d.py

def __init__(self, sampler):
    r"""
    Args:
        sampler (Union[AbstractDiscreteDistribution,AbstractTrueMeasure]): Either  

            - a discrete distribution from which to transform samples, or
            - a true measure by which to compose a transform.
    """
    self.sampler = sampler 
    assert self.sampler.d==2
    self.true_measure = Uniform(self.sampler,lower_bound=[-4,-3],upper_bound=[7,8])
    super(Multimodal2d,self).__init__(dimension_indv=(),dimension_comb=(),parallel=False)

Hartmann6d

Bases: AbstractIntegrand

Wrapper around BoTorch's implementation of the Augmented Hartmann function in dimension \(d=6\).

Examples:

>>> integrand = Hartmann6d(DigitalNetB2(6,seed=7))
>>> y = integrand(2**10)
>>> print("%.4f"%y.mean())
-0.2644
>>> integrand.true_measure
Uniform (AbstractTrueMeasure)
    lower_bound     0
    upper_bound     1

With independent replications

>>> integrand = Hartmann6d(DigitalNetB2(6,seed=7,replications=2**4))
>>> y = integrand(2**6)
>>> y.shape
(16, 64)
>>> muhats = y.mean(-1) 
>>> muhats.shape 
(16,)
>>> print("%.4f"%muhats.mean())
-0.2599

Parameters:

Name	Type	Description	Default
sampler	Union[AbstractDiscreteDistribution, AbstractTrueMeasure]	Either a discrete distribution from which to transform samples, or a true measure by which to compose a transform.	required

Source code in qmcpy/integrand/hartmann6d.py

def __init__(self, sampler):
    r"""
    Args:
        sampler (Union[AbstractDiscreteDistribution,AbstractTrueMeasure]): Either  

            - a discrete distribution from which to transform samples, or
            - a true measure by which to compose a transform.
    """
    self.sampler = sampler
    assert self.sampler.d==6
    self.true_measure = Uniform(self.sampler,lower_bound=0,upper_bound=1)
    super(Hartmann6d,self).__init__(dimension_indv=(),dimension_comb=(),parallel=False)
    from botorch.test_functions.multi_fidelity import AugmentedHartmann
    self.ah = AugmentedHartmann(negate=False)

FourBranch2d

Bases: AbstractIntegrand

Four Branch function in \(d=2\).

\[g(\boldsymbol{t}) = \min \begin{cases} 3+0.1(t_0-t_1)^2-\frac{t_0-t_1}{\sqrt{2}} \\ 3+0.1(t_0-t_1)^2+\frac{t_0-t_1}{\sqrt{2}} \\ t_0-t_1 + 7/\sqrt{2} \\ t_1-t_0 + 7/\sqrt{2}\end{cases}, \qquad \boldsymbol{T}=(T_0,T_1) \sim \mathcal{U}[-8,8]^2.\]

Examples:

>>> integrand = FourBranch2d(DigitalNetB2(2,seed=7))
>>> y = integrand(2**10)
>>> print("%.4f"%y.mean())
-2.4995
>>> integrand.true_measure
Uniform (AbstractTrueMeasure)
    lower_bound     -8
    upper_bound     2^(3)

With independent replications

>>> integrand = FourBranch2d(DigitalNetB2(2,seed=7,replications=2**4))
>>> y = integrand(2**6)
>>> y.shape
(16, 64)
>>> muhats = y.mean(-1) 
>>> muhats.shape 
(16,)
>>> print("%.4f"%muhats.mean())
-2.5042

Parameters:

Name	Type	Description	Default
sampler	Union[AbstractDiscreteDistribution, AbstractTrueMeasure]	Either a discrete distribution from which to transform samples, or a true measure by which to compose a transform.	required

Source code in qmcpy/integrand/fourbranch2d.py

def __init__(self, sampler):
    r"""
    Args:
        sampler (Union[AbstractDiscreteDistribution,AbstractTrueMeasure]): Either  

            - a discrete distribution from which to transform samples, or
            - a true measure by which to compose a transform.
    """
    self.sampler = sampler
    assert self.sampler.d==2
    self.true_measure = Uniform(self.sampler,lower_bound=-8,upper_bound=8)
    super(FourBranch2d,self).__init__(dimension_indv=(),dimension_comb=(),parallel=False)

BayesianLRCoeffs

Bases: AbstractIntegrand

Logistic Regression Coefficients computed as the posterior mean in a Bayesian framework.

Examples:

>>> integrand = BayesianLRCoeffs(DigitalNetB2(3,seed=7),feature_array=np.arange(8).reshape((4,2)),response_vector=[0,0,1,1])
>>> y = integrand(2**10)
>>> y.shape
(2, 3, 1024)
>>> y.mean(-1)
array([[ 0.04517466, -0.01103669, -0.06614381],
       [ 0.02162049,  0.02162049,  0.02162049]])

With independent replications

>>> integrand = BayesianLRCoeffs(DigitalNetB2(3,seed=7,replications=2**4),feature_array=np.arange(8).reshape((4,2)),response_vector=[0,0,1,1])
>>> y = integrand(2**6)
>>> y.shape
(2, 3, 16, 64)
>>> muhats = y.mean(-1) 
>>> muhats.shape 
(2, 3, 16)
>>> muhats.mean(-1)
array([[ 0.0587368 , -0.01718134, -0.07203021],
       [ 0.02498059,  0.02498059,  0.02498059]])

Parameters:

Name	Type	Description	Default
sampler	Union[AbstractDiscreteDistribution, AbstractTrueMeasure]	Either a discrete distribution from which to transform samples, or a true measure by which to compose a transform.	required
feature_array	ndarray	Array of features with shape \((N,d-1)\) where \(N\) is the number of observations and \(d\) is the dimension.	required
response_vector	ndarray	Binary responses vector of length \(N\).	required
prior_mean	ndarray	Length \(d\) vector of prior means, one for each coefficient. The first \(d-1\) inputs correspond to the \(d-1\) features. The last input corresponds to the intercept coefficient.	0
prior_covariance	ndarray	Prior covariance array with shape \((d,d)\) d x d where indexing is consistent with the prior mean.	10

Source code in qmcpy/integrand/bayesian_lr_coeffs.py

def __init__(self, sampler, feature_array, response_vector, prior_mean=0, prior_covariance=10):
    r"""
    Args:
        sampler (Union[AbstractDiscreteDistribution,AbstractTrueMeasure]): Either  

            - a discrete distribution from which to transform samples, or
            - a true measure by which to compose a transform.
        feature_array (np.ndarray): Array of features with shape $(N,d-1)$ where $N$ is the number of observations and $d$ is the dimension.
        response_vector (np.ndarray): Binary responses vector of length $N$.
        prior_mean (np.ndarray): Length $d$ vector of prior means, one for each coefficient.

            - The first $d-1$ inputs correspond to the $d-1$ features. 
            - The last input corresponds to the intercept coefficient.
        prior_covariance (np.ndarray): Prior covariance array with shape $(d,d)$ d x d where indexing is consistent with the prior mean.
    """
    self.prior_mean = prior_mean
    self.prior_covariance = prior_covariance
    self.sampler = sampler
    self.true_measure = Gaussian(self.sampler, mean=self.prior_mean, covariance=self.prior_covariance)
    self.feature_array = np.array(feature_array,dtype=float)
    self.response_vector = np.array(response_vector,dtype=float)
    obs,dm1 = self.feature_array.shape
    self.num_coeffs = dm1+1
    if self.num_coeffs!=self.true_measure.d:
        ParameterError("sampler must have dimension one more than the number of features in the feature_array.")
    if self.response_vector.shape!=(obs,) or ((self.response_vector!=0)&(self.response_vector!=1)).any():
        ParameterError("response_vector must have the same length as feature_array and contain only 0 or 1 entries.")
    self.feature_array = np.column_stack((self.feature_array,np.ones((obs,1))))
    super(BayesianLRCoeffs,self).__init__(dimension_indv=(2,self.num_coeffs),dimension_comb=self.num_coeffs,parallel=False)

Sin1d

Bases: AbstractIntegrand

Sine function in \(d=1\) dimension.

\[g(t) = \sin(t), \qquad t \sim \mathcal{U}[0,2\pi k]\]

Examples:

>>> integrand = Sin1d(DigitalNetB2(1,seed=7),k=1)
>>> y = integrand(2**10)
>>> print("%.4e"%y.mean())
-1.3582e-10
>>> integrand.true_measure
Uniform (AbstractTrueMeasure)
    lower_bound     0
    upper_bound     6.283

With independent replications

>>> integrand = Sin1d(DigitalNetB2(1,seed=7,replications=2**4),k=1)
>>> y = integrand(2**6)
>>> y.shape
(16, 64)
>>> muhats = y.mean(-1) 
>>> muhats.shape 
(16,)
>>> print("%.4e"%muhats.mean())
7.0800e-04

Parameters:

Name	Type	Description	Default
sampler	Union[AbstractDiscreteDistribution, AbstractTrueMeasure]	Either a discrete distribution from which to transform samples, or a true measure by which to compose a transform.	required
k	float	The true measure will be uniform between \(0\) and \(2 \pi k\).	1

Source code in qmcpy/integrand/sin1d.py

def __init__(self, sampler, k=1):
    r"""
    Args:
        sampler (Union[AbstractDiscreteDistribution,AbstractTrueMeasure]): Either  

            - a discrete distribution from which to transform samples, or
            - a true measure by which to compose a transform.
        k (float): The true measure will be uniform between $0$ and $2 \pi k$.
    """
    self.sampler = sampler
    self.k = k
    assert self.sampler.d==1
    self.true_measure = Uniform(self.sampler,lower_bound=0,upper_bound=2*self.k*np.pi)
    super(Sin1d,self).__init__(dimension_indv=(),dimension_comb=(),parallel=False)

Linear0

Bases: AbstractIntegrand

Linear Function with analytic mean \(0\).

\[g(\boldsymbol{t}) = \sum_{j=1}^d t_j \qquad \boldsymbol{T} \sim \mathcal{U}[0,1]^d.\]

Examples:

>>> integrand = Linear0(DigitalNetB2(100,seed=7))
>>> y = integrand(2**10)
>>> print("%.4e"%y.mean())
6.0560e-05

With independent replications

>>> integrand = Linear0(DigitalNetB2(100,seed=7,replications=2**4))
>>> y = integrand(2**6)
>>> y.shape
(16, 64)
>>> muhats = y.mean(-1) 
>>> muhats.shape 
(16,)
>>> print("%.4e"%muhats.mean())
-9.8203e-05

Parameters:

Name	Type	Description	Default
sampler	Union[AbstractDiscreteDistribution, AbstractTrueMeasure]	Either a discrete distribution from which to transform samples, or a true measure by which to compose a transform.	required

Source code in qmcpy/integrand/linear0.py

def __init__(self, sampler):
    r"""
    Args:
        sampler (Union[AbstractDiscreteDistribution,AbstractTrueMeasure]): Either  

            - a discrete distribution from which to transform samples, or
            - a true measure by which to compose a transform.
    """
    self.sampler = sampler
    self.true_measure = Uniform(self.sampler, lower_bound=-.5, upper_bound=.5)
    super(Linear0,self).__init__(dimension_indv=(),dimension_comb=(),parallel=False)

UML Specific