GAIL---Guaranteed Automatic Integration Library in MATLAB: Documentation for Version 2.1

Automatic and adaptive approximation, optimization, or integration of functions in a cone with guarantee of accuracy is a relatively new paradigm. Our purpose is to create an open-source MATLAB package, Guaranteed Automatic Integration Library (GAIL), following the philosophy of reproducible research and sustainable practices of robust scientific software development. For our conviction that true scholarship in computational sciences are characterized by reliable reproducibility, we employ the best practices in mathematical research and software engineering known to us and available in MATLAB. This document describes the key features of functions in GAIL, which includes one-dimensional function approximation and minimization using linear splines, one-dimensional numerical integration using trapezoidal rule, and last but not least, mean estimation and multidimensional integration by Monte Carlo methods or Quasi Monte Carlo methods.


Introduction
Automatic and adaptive approximation, optimization, or integration of functions in a cone with guarantee of accuracy is a relatively new paradigm [7]. Our purpose is to create an open-source MATLAB package, Guaranteed Automatic Integration Library (GAIL) [5], following the philosophy of reproducible research championed by Claerbout [6] and Donoho [1], and sustainable practices of robust scientific software development [12]. For our conviction that true scholarship in computational sciences are characterized by reliable reproducibility [3,4,2], we employ the best practices in mathematical research and software engineering known to us and available in MATLAB.
The rest of this document describes the key features of functions in GAIL, which includes one-dimensional 1.5 Installation Instruction 1. Unzip the contents of the zip file to a directory and maintain the existing directory and subdirectory structure. (Please note: If you install into the toolbox subdirectory of the MATLAB program hierarchy, you will need to click the button "Update toolbox path cache" from the File/Preferences... dialog in MATLAB.) 2. In MATLAB, add the GAIL directory to your path. This can be done by running GAIL Install.m. Alternatively, this can be done by selecting "File/Set Path..." from the main or Command window menus, or with the command pathtool. We recommend that you select the "Save" button on this dialog so that GAIL is on the path automatically in future MATLAB sessions. 3. To check if you have installed GAIL successfully, type help funappx g to see if its documentation shows up.
Alternatively, you could do this: 1. Download DownloadInstallGail 2 1.m and put it where you want GAIL to be installed.

Execute it in MATLAB.
To uninstall GAIL, execute GAIL Uninstall.
To reinstall GAIL, execute GAIL Install.

Tests
We provide quick doctests for each of the functions above. To run doctests in funappx g, for example, issue the command doctest funappx g.
We also provide unit tests for MATLAB version 8 or later. To run unit tests for funmin g, for instance, execute run(ut funmin g).

Contact Information
Please send any queries, questions, or comments to gail-users@googlegroups.com

Website
For more information about GAIL, visit GAIL Project website.

Description
fappx = funappx g(f) approximates function f on the default interval [0,1] by an approximated function handle fappx within the guaranteed absolute error tolerance of 1e-6. When Matlab version is higher or equal to 8.3, fappx is an interpolant generated by griddedInterpolant. When Matlab version is lower than 8.3, fappx is a function handle generated by ppval and interp1. Input f is a function handle. The statement y = f(x) should accept a vector argument x and return a vector y of function values that is of the same size as x. fappx = funappx g(f,a,b,abstol) for a given function f and the ordered input parameters that define the finite interval [a,b], and a guaranteed absolute error tolerance abstol. fappx = funappx g(f,'a',a,'b',b,'abstol',abstol) approximates function f on the finite interval [a,b], given a guaranteed absolute error tolerance abstol. All four field-value pairs are optional and can be supplied in different order. fappx = funappx g(f,in param) approximates function f on the finite interval [in param.a,in param.b], given a guaranteed absolute error tolerance in param.abstol. If a field is not specified, the default value is used.
Input Arguments If the function to be approximated, f satisfies the cone condition for each sub interval [t l−1 , t l ], where 1 ≤ l ≤ L, then the f appx |output by this algorithm is guaranteed to satisfy f − f appx ∞ ≤ abstol.  [fmin, out param] = funmin g(f,...) returns minimum value fmin of function f and an output structure out param.

Optional Input Arguments
• in param.nlo -lower bound of initial number of points we used, default value is 10 • in param.nhi -upper bound of initial number of points we used, default value is 1000 • in param.nmax -cost budget, default value is 1e7.

Output Arguments
• out param.f -input function

Guarantee
If the function to be minimized, f satisfies the cone condition then the fmin output by this algorithm is guaranteed to satisfy provided the flag exitflag = 0.   [q, out param] = integral g(f,...) returns the approximated integration q and output structure out param.

Description
Input Arguments • out param.tauchange -it is true if the cone constant has been changed, false otherwise. See [1] for details. If true, you may wish to change the input in param.ninit to a larger number.
• out param.tauchange -it is true if the cone constant has been changed, false otherwise. See [1] for details. If true, you may wish to change the input in param.ninit to a larger number.
• out param.iter -number of iterations • out param.npoints -number of points we need to reach the guaranteed absolute error tolerance abstol.
• out param.errest -approximation error defined as the differences between the true value and the approximated value of the integral.
• out param.nstar -final value of the parameter defining the cone of functions for which this algorithm is guaranteed; nstar = ninit-2 initially and is increased as necessary

Guarantee
If the function to be integrated, f satisfies the cone condition , then the q output by this algorithm is guaranteed to satisfy

Description
tmu = meanMC g(Yrand) estimates the mean, mu, of a random variable Y to within a specified generalized error tolerance, tolfun:=max(abstol,reltol*| mu |), i.e., | mu -tmu | <= tolfun with probability at least 1alpha, where abstol is the absolute error tolerance, and reltol is the relative error tolerance. Usually the reltol determines the accuracy of the estimation, however, if the | mu | is rather small, the abstol determines the accuracy of the estimation. The default values are abstol=1e-2, reltol=1e-1, and alpha=1%. Input Yrand is a function handle that accepts a positive integer input n and returns an n x 1 vector of IID instances of the random variable Y. tmu = meanMC g(Yrand,abstol,reltol,alpha) estimates the mean of a random variable Y to within a specified generalized error tolerance tolfun with guaranteed confidence level 1-alpha using all ordered parsing inputs abstol, reltol, alpha. tmu = meanMC g(Yrand,'abstol',abstol,'reltol',reltol,'alpha',alpha) estimates the mean of a random variable Y to within a specified generalized error tolerance tolfun with guaranteed confidence level 1-alpha. All the field-value pairs are optional and can be supplied in different order, if a field is not supplied, the default value is used.
[tmu, out param] = meanMC g(Yrand,in param) estimates the mean of a random variable Y to within a specified generalized error tolerance tolfun with the given parameters in param and produce the estimated mean tmu and output parameters out param. If a field is not specified, the default value is used.

Input Arguments
• Yrand -the function for generating n IID instances of a random variable Y whose mean we want to estimate. Y is often defined as a function of some random variable X with a simple distribution. The input of Yrand should be the number of random variables n, the output of Yrand should be n function values. For example, if Y = X.ˆ2 where X is a standard uniform random variable, then one may define Yrand = @(n) rand(n,1).ˆ2.
• in param.abstol -the absolute error tolerance, which should be positive, default value is 1e-2.
• in param.reltol -the relative error tolerance, which should be between 0 and 1, default value is 1e-1.
• in param.alpha -the uncertainty, which should be a small positive percentage. default value is 1%.

Optional Input Arguments
• in param.fudge -standard deviation inflation factor, which should be larger than 1, default value is 1.2.
• in param.nSig -initial sample size for estimating the sample variance, which should be a moderate large integer at least 30, the default value is 1e4.
• in param.n1 -initial sample size for estimating the sample mean, which should be a moderate large positive integer at least 30, the default value is 1e4.
• in param.tbudget -the time budget in seconds to do the two-stage estimation, which should be positive, the default value is 100 seconds.
• in param.nbudget -the sample budget to do the two-stage estimation, which should be a large positive integer, the default value is 1e9.

Output Arguments
• tmu -the estimated mean of Y.
• out param.tau -the iteration step.
• out param.n -the sample size used in each iteration.
• out param.nremain -the remaining sample budget to estimate mu. It was calculated by the sample left and time left.
• out param.hmu -estimated mean in each iteration.
• out param.tol -the reliable upper bound on error for each iteration.
• out param.var -the sample variance.
• out param.exit -the state of program when exiting. • out param.time -the time elapsed in seconds.
• out param.flag -parameter checking status 1 checked by meanMC g

Guarantee
This algorithm attempts to calculate the mean, mu, of a random variable to a prescribed error tolerance, tolfun:= max(abstol,reltol*| mu |), with guaranteed confidence level 1-alpha. If the algorithm terminated without showing any warning messages and provide an answer tmu, then the follow inequality would be satisfied: Pr(| mu -tmu | <= tolfun) >= 1-alpha The cost of the algorithm, N tot, is also bounded above by N up, which is defined in terms of abstol, reltol, nSig, n1, fudge, kurtmax, beta. And the following inequality holds: Pr (N tot <= N up) >= 1-beta Please refer to our paper for detailed arguments and proofs.

meanMCBer g
Monte Carlo method to estimate the mean of a Bernoulli random variable to within a specified absolute error tolerance with guaranteed confidence level 1-alpha.

Description
pHat = meanMCBer g(Yrand) estimates the mean of a Bernoulli random variable Y to within a specified absolute error tolerance with guaranteed confidence level 99%. Input Yrand is a function handle that accepts a positive integer input n and returns a n x 1 vector of IID instances of the Bernoulli random variable Y. pHat = meanMCBer g(Yrand,abstol,alpha,nmax) estimates the mean of a Bernoulli random variable Y to within a specified absolute error tolerance with guaranteed confidence level 1-alpha using all ordered parsing inputs abstol, alpha and nmax. pHat = meanMCBer g(Yrand,'abstol',abstol,'alpha',alpha,'nmax',nmax) estimates the mean of a Bernoulli random variable Y to within a specified absolute error tolerance with guaranteed confidence level 1-alpha. All the field-value pairs are optional and can be supplied in different order.
[pHat, out param] = meanMCBer g(Yrand,in param) estimates the mean of a Bernoulli random variable Y to within a specified absolute error tolerance with the given parameters in param and produce the estimated mean pHat and output parameters out param.

Input Arguments
• Yrand -the function for generating IID instances of a Bernoulli random variable Y whose mean we want to estimate.
• pHat -the estimated mean of Y.
• in param.abstol -the absolute error tolerance, the default value is 1e-2.
• in param.alpha -the uncertainty, the default value is 1%.
• in param.nmax -the sample budget, the default value is 1e9.

Output Arguments
• out param.n -the total sample used.
• out param.time -the time elapsed in seconds.
• out param.exit -the state of program when exiting.
0 Success 1 Not enough samples to estimate p with guarantee

Guarantee
If the sample size is calculated according Hoeffding's inequality, which equals to ceil(log(2/out param.alpha)/(2*out param.abstolˆ2)), then the following inequality must be satisfied: Pr(| p -pHat | <= abstol) >= 1-alpha. Here p is the true mean of Yrand, and pHat is the output of MEANMCBER G. Also, the cost is deterministic.

Description
[Q,out param] = cubMC g(f,hyperbox) estimates the integral of f over hyperbox to within a specified generalized error tolerance, tolfun = max(abstol, reltol*| I |), i.e., | I -Q | <= tolfun with probability at least 1-alpha, where abstol is the absolute error tolerance, and reltol is the relative error tolerance. Usually the reltol determines the accuracy of the estimation, however, if the | I | is rather small, the abstol determines the accuracy of the estimation. The default values are abstol=1e-2, reltol=1e-1, and alpha=1%. Input f is a function handle that accepts an n x d matrix input, where d is the dimension of the hyperbox, and n is the number of points being evaluated simultaneously. The input hyperbox is a 2 x d matrix, where the first row corresponds to the lower limits and the second row corresponds to the upper limits. Q = cubMC g(f,hyperbox,measure,abstol,reltol,alpha) estimates the integral of function f over hyperbox to within a specified generalized error tolerance tolfun with guaranteed confidence level 1-alpha using all ordered parsing inputs f, hyperbox, measure, abstol, reltol, alpha, fudge, nSig, n1, tbudget, nbudget, flag. The input f and hyperbox are required and others are optional. Q = cubMC g(f,hyperbox,'measure',measure,'abstol',abstol,'reltol',reltol,'alpha',alpha) estimates the integral of f over hyperbox to within a specified generalized error tolerance tolfun with guaranteed confidence level 1-alpha. All the field-value pairs are optional and can be supplied in different order. If an input is not specified, the default value is used.
[Q out param] = cubMC g(f,hyperbox,in param) estimates the integral of f over hyperbox to within a specified generalized error tolerance tolfun with the given parameters in param and produce output parameters out param and the integral Q.

Input Arguments
• f -the integrand.
• in param.measure -the measure for generating the random variable, the default is 'uniform'. The other measure could be handled is 'normal'/'Gaussian'. The input should be a string type, hence with quotes.
• in param.abstol -the absolute error tolerance, the default value is 1e-2.
• in param.reltol -the relative error tolerance, the default value is 1e-1.
• in param.alpha -the uncertainty, the default value is 1%.

Optional Input Arguments
• in param.fudge -the standard deviation inflation factor, the default value is 1.2.
• in param.nSig -initial sample size for estimating the sample variance, which should be a moderate large integer at least 30, the default value is 1e4.
• in param.n1 -initial sample size for estimating the sample mean, which should be a moderate large positive integer at least 30, the default value is 1e4.
• in param.tbudget -the time budget to do the estimation, the default value is 100 seconds.
• in param.nbudget -the sample budget to do the estimation, the default value is 1e9.
• in param.flag -the value corresponds to parameter checking status.
0 not checked 1 checked by meanMC g 2 checked by cubMC g

Output Arguments
• Q -the estimated value of the integral.
• out param.n -the sample size used in each iteration.
• out param.nremain -the remaining sample budget to estimate I. It was calculated by the sample left and time left.
• out param.tau -the iteration step.
• out param.hmu -estimated integral in each iteration.
• out param.tol -the reliable upper bound on error for each iteration.
• out param.kurtmax -the upper bound on modified kurtosis.
• out param.time -the time elapsed in seconds.
• out param.var -the sample variance.
• out param.exit -the state of program when exiting.

Guarantee
This algorithm attempts to calculate the integral of function f over a hyperbox to a prescribed error tolerance tolfun:= max(abstol,reltol*| I |) with guaranteed confidence level 1-alpha. If the algorithm terminated without showing any warning messages and provide an answer Q, then the follow inequality would be satisfied: The cost of the algorithm, N tot, is also bounded above by N up, which is a function in terms of abstol, reltol, nSig, n1, fudge, kurtmax, beta. And the following inequality holds: Pr (N tot <= N up) >= 1-beta Please refer to our paper for detailed arguments and proofs.

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

Description
[q,out param] = cubLattice g(f,hyperbox) estimates the integral of f over the d-dimensional region described by hyperbox, and with an error guaranteed not to be greater than a specific generalized error tolerance, tolfun:=max(abstol,reltol*| integral(f) |). Input f is a function handle. f should accept an n x d matrix input, where d is the dimension and n is the number of points being evaluated simultaneously. The input hyperbox is a 2 x d matrix, where the first row corresponds to the lower limits and the second row corresponds to the upper limits of the integral. Given the construction of our Lattices, d must be a positive integer with 1<=d<=250. q = cubLattice g(f,hyperbox,measure,abstol,reltol) estimates the integral of f over the hyperbox. The answer is given within the generalized error tolerance tolfun. All parameters should be input in the order specified above. If an input is not specified, the default value is used. Note that if an input is not specified, the remaining tail cannot be specified either. Inputs f and hyperbox are required. The other optional inputs are in the correct order: measure,abstol,reltol,shift,mmin,mmax,fudge,transform,toltype and theta. q = cubLattice g(f,hyperbox,'measure',measure,'abstol',abstol,'reltol',reltol) estimates the integral of f over the hyperbox. The answer is given within the generalized error tolerance tolfun. All the field-value pairs are optional and can be supplied in any order. If an input is not specified, the default value is used. q = cubLattice g(f,hyperbox,in param) estimates the integral of f over the hyperbox. The answer is given within the generalized error tolerance tolfun.

Input Arguments
• f -the integrand whose input should be a matrix n x d where n is the number of data points and d the dimension, which cannot be greater than 250. By default f is f=@ x.ˆ2.
• hyperbox -the integration region defined by its bounds. It must be a 2 x d matrix, where the first row corresponds to the lower limits and the second row corresponds to the upper limits of the integral. The default value is [0;1].
• in param.measure -for f(x)*mu(dx), we can define mu(dx) to be the measure of a uniformly distributed random variable in they hyperbox or normally distributed with covariance matrix I d. The only possible values are 'uniform' or 'normal'. For 'uniform', the hyperbox must be a finite volume while for 'normal', the hyperbox can only be defined as (-Inf,Inf)ˆd. By default it is 'uniform'.

Optional Input Arguments
• in param.shift -the Rank-1 lattices can be shifted to avoid the origin or other particular points. By default we consider a uniformly [0,1) random shift.
• in param.mmin -the minimum number of points to start is 2ˆmmin. The cone condition on the Fourier coefficients decay requires a minimum number of points to start. The advice is to consider at least mmin=10. mmin needs to be a positive integer with mmin<=mmax. By default it is 10.
• in param.mmax -the maximum budget is 2ˆmmax. By construction of our Lattices generator, mmax is a positive integer such that mmin<=mmax<=26. The default value is 24.
• in param.fudge -the positive function multiplying the finite sum of Fast Fourier coefficients specified in the cone of functions. This input is a function handle. The fudge should accept an array of nonnegative integers being evaluated simultaneously. For more technical information about this parameter, refer to the references. By default it is @(m) 5*2.ˆ-m. 'C1sin' : Sidi's transform with sine, preserving the first derivative. This is in general a better option than 'C1'.
• in param.toltype -this is the generalized tolerance function. There are two choices, 'max' which takes max(abstol,reltol*| integral(f) | ) and 'comb' which is the linear combination theta*abstol+(1theta)*reltol*| integral(f) | . Theta is another parameter to be specified with 'comb'(see below). For pure absolute error, either choose 'max' and set reltol = 0 or choose 'comb' and set theta = 1. For pure relative error, either choose 'max' and set abstol = 0 or choose 'comb' and set theta = 0. Note that with 'max', the user can not input abstol = reltol = 0 and with 'comb', if theta = 1 abstol con not be 0 while if theta = 0, reltol can not be 0. By default toltype is 'max'.
• in param.theta -this input is parametrizing the toltype 'comb'. Thus, it is only active when the toltype chosen is 'comb'. It establishes the linear combination weight between the absolute and relative tolerances theta*abstol+(1-theta)*reltol*| integral(f) |. Note that for theta = 1, we have pure absolute tolerance while for theta = 0, we have pure relative tolerance. By default, theta=1.

Output Arguments
• q -the estimated value of the integral.
• out param.d -dimension over which the algorithm integrated.
• out param.n -number of Rank-1 lattice points used for computing the integral of f.
• out param.bound err -predicted bound on the error based on the cone condition. If the function lies in the cone, the real error will be smaller than generalized tolerance.
• out param.time -time elapsed in seconds when calling cubLattice g.

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

Description
[q,out param] = cubSobol g(f,hyperbox) estimates the integral of f over the d-dimensional region described by hyperbox, and with an error guaranteed not to be greater than a specific generalized error tolerance, tolfun:=max(abstol,reltol*| integral(f) |). Input f is a function handle. f should accept an n x d matrix input, where d is the dimension and n is the number of points being evaluated simultaneously. The input hyperbox is a 2 x d matrix, where the first row corresponds to the lower limits and the second row corresponds to the upper limits of the integral. Given the construction of Sobol' sequences, d must be a positive integer with 1<=d<=1111. q = cubSobol g(f,hyperbox,measure,abstol,reltol) estimates the integral of f over the hyperbox. The answer is given within the generalized error tolerance tolfun. All parameters should be input in the order specified above. If an input is not specified, the default value is used. Note that if an input is not specified, the remaining tail cannot be specified either. Inputs f and hyperbox are required. The other optional inputs are in the correct order: measure,abstol,reltol,mmin,mmax,fudge,toltype and theta. q = cubSobol g(f,hyperbox,'measure',measure,'abstol',abstol,'reltol',reltol) estimates the integral of f over the hyperbox. The answer is given within the generalized error tolerance tolfun. All the field-value pairs are optional and can be supplied in any order. If an input is not specified, the default value is used. q = cubSobol g(f,hyperbox,in param) estimates the integral of f over the hyperbox. The answer is given within the generalized error tolerance tolfun.

Input Arguments
• f -the integrand whose input should be a matrix n x d where n is the number of data points and d the dimension, which cannot be greater than 1111. By default f is f=@ x.ˆ2.
• hyperbox -the integration region defined by its bounds. It must be a 2 x d matrix, where the first row corresponds to the lower limits and the second row corresponds to the upper limits of the integral. The default value is [0;1].
• in param.measure -for f(x)*mu(dx), we can define mu(dx) to be the measure of a uniformly distributed random variable in the hyperbox or normally distributed with covariance matrix I d. The only possible values are 'uniform' or 'normal'. For 'uniform', the hyperbox must be a finite volume while for 'normal', the hyperbox can only be defined as (-Inf,Inf)ˆd. By default it is 'uniform'.