We will use the package to test Monte Carlo integration over the unit hypercube \(C_n=[0,1]^n\). The package covers test functions for five different domains, but the interfaces are the same in every case. For each function in this package the exact value of the integral over the respective integration domain is known for any dimension \(n\geq 1\). This allows to check and contrast multivariate integration routines for in different settings such as low, medium and high dimensional. In this vignette we use as an example Monte Carlo type integration due to its simplicity. The interface to the package is the same in all cases:
We need to instantiate an object from one of the classes. This object represents the function that we want to integrate. Each object (function) has a mandatory variable that represents the dimension of the integral. Other parameters depend on the specific functions. With the instance we then have access to basic information about the function and can evaluate it for a given set of points. These evaluation points are arranged in a matrix where each row represents an evaluation point. The number of columns needs to be equal to the dimension variable that is set when creating the object (function). To check which points in such a matrix are in the integration domain linked to the function one can use checkDomain. In total there are six methods that can be used for an instance.
Warning: The numerical integration scheme utilized in the following section is only used to showcase the package and its functions. It is not intended as an examples of a “best practice” when it comes to multivariate numerical integration.
We exemplify the usage of the package by numerically integrate over the unit hypercube \(C_n = [0,1]^n\). In this case Monte Carlo integration works straight out of the box with standard uniform (pseudo) random numbers. If \(\vec{U} = (U_1,\ldots,U_n)\) is uniformly distributed over \(C_n\) we can approximate the integral of a function \(f\) via \[ \int_{C_n} f(\vec{x}) d\vec{x} = \mathbb{E}[f(\vec{U})] \approx \frac{1}{n} \sum_{i=1}^n f(\vec{U}_i), \] where the law of large numbers is used with independent copies \(\vec{U}_i\), \(1\leq i\leq n\), of \(\vec{U}\). It is easy to generate such independent copies by simply drawing independent (pseudo) random numbers. For example, for the specific choice of \(f\) as \[ f(\vec{x}) = (\cos( a_1x_1+\ldots+a_n x_n ))^2 \] we have \[ \int_{C_n} f(\vec{x}) d\vec{x} = \frac{1}{2} + \frac{1}{2}\cos(a_1+\ldots+a_n)\prod_{i=1}^n\frac{\sin(a_j)}{a_j}, \] where \(a_i\neq 0\). We can now compare the numerical integration outlined above to the theoretical value as follows.
require(multIntTestFunc)
set.seed(42)
n <- as.integer(4) #the dimension needs to be of type "integer"
coeffs <- 1:n #get example coefficients
f <- new("unitCube_cos2",coeffs=coeffs,dim=n) #create instance
#print some information on the function
print(getIntegrationDomain(f))## [1] "standard unit hypercube: 0 <= x_i <= 1"print(getTags(f))## [1] "unit hypercube" "smooth" "cos" "continuous"print(getReferences(f))## [1] "C.1"#evaluate integral by simple Monte Carlo integration
N <- 10000 #choose number of node points
U <- matrix(runif(N*n),N,n) #arrange points in matrix
D <- domainCheck(f,U) #check if the points are in the integration domain
U <- U[D,] #select the points in the integration domain (here f is supported over C_n so we take all points)
U <- as.matrix(U,ncol=n) #arrange points in matrix - only necessary because for n=1 we might get a vector, but we need a matrix input for our evaluate function
eval <- evaluate(f,U) #evaluate the function f on the points U
approx <- mean(eval) #approximate the expectation by average
print(approx)## [1] 0.4967258#print exact integral
print(exactIntegral(f))## [1] 0.5014285