Monte-Carlo and Quasi Monte-Carlo Integration

Execute Monte-Carlo or Quasi Monte-Carlo Integration for given integrand function.

There are two pre-defined Low discrepancy point sets, Niederreiter-Xing point set and Sobol point set. Both of them are selected to have low WAFOM.

Integrand function should recieve numeric vector and shoud return numeic value. Points passed to integrand are in the range (0, 1)^s. Returned value of integration is mean of returned value of integrand, not summation. To get summation, multiply by 2^m.

First, define an example integrand function; calculate n-sphere (0,1)^n.

library(rmcqmcint)
unit.nsphere <- function(point) {
    if (sqrt(sum(point^2)) <= 1.0) {
        r <- 1.0
    } else {
        r <- 0.0
    }
    return(r)
}

example 1: Calculate 4d-sphare by qmc, with default digital net.

rs <- qmcint(integrand=unit.nsphere, N=100, s=4)
rs$mean

## [1] 0.3100977

rs$absError

## [1] 0.006666557

1 / 2 * pi^2 / 16 # to compare true value.

## [1] 0.3084251

example 2: Calculate 5d-sphare by qmc, with sobol sequence low WAFOM digital net.

rs <- qmcint(integrand=unit.nsphere, N=100, s=5, digitalNetID = 2)
rs$mean

## [1] 0.164375

rs$absError

## [1] 0.001514922

8 / 15 * pi^2 / 32 # to compare true value.

## [1] 0.1644934

example 3: Calculate 4d-sphare by mc.

rs <- qmcint(integrand=unit.nsphere, N=100, s=4)
rs$mean

## [1] 0.3100977

rs$absError

## [1] 0.006666557

1 / 2 * pi^2 / 16 # to compare true value.

## [1] 0.3084251

example 4: Calculate 5d-sphare by mc.

rs <- qmcint(integrand=unit.nsphere, N=100, s=5)
rs$mean

## [1] 0.1638379

rs$absError

## [1] 0.002203611

8 / 15 * pi^2 / 32 # to compare true value.

## [1] 0.1644934

Define another example integrand function; Oscillatory function in Genz test functions.

make.oscillatory <- function(a) {
    oscillatory <- function(point) {
        return(cos(sum(a * point)))
    }
    return(oscillatory)
}

If sum(a) is pi, integration of this function in (0, 1)^s will become zero. ID=3 means to use Sobol sequence up to dimention 21201.

n <- 100
id <- 3
s <- 1000
m <- 10
p <- 0.99
a <- rep(pi / s, length=s)
osc <- make.oscillatory(a)
rs <- qmcint(osc, N=n, s=s, digitalNetID=id, m=m, probability=p)
rs$mean

## [1] 1.167836e-06

rs$absError

## [1] 0.0002671781