sphunif (version 1.4.3)

int_sph_MC: Monte Carlo integration of functions on the (hyper)sphere

Description

Monte Carlo approximation of the integral $$\int_{S^{p-1}}f(x)\,\mathrm{d}x$$ of a function \(f:S^{p-1} \rightarrow R\) defined on the (hyper)sphere \(S^{p-1}:=\{{\bf x}\in R^p:||{\bf x}||=1\}\), \(p\ge 2\).

Usage

int_sph_MC(f, p, M = 10000, cores = 1, chunks = ceiling(M/1000),
  seeds = NULL, ...)

Value

A scalar with the approximate integral.

Arguments

f

function to be integrated. Its first argument must be the (hyper)sphere position. Must be vectorized and return a vector of size nrow(x) for a matrix input x. See examples.

p

integer giving the dimension of the ambient space \(R^p\) that contains \(S^{p-1}\).

M

number of Monte Carlo samples. Defaults to 1e4.

cores

number of cores to perform the integration. Defaults to 1.

chunks

number of chunks to split the M Monte Carlo samples. Useful for parallelizing the integration in chunks tasks containing ceiling(M / chunks) replications. Useful also for avoiding memory bottlenecks when M is large. Defaults to
ceiling(M / 1e3).

seeds

if provided, a vector of size chunks for fixing the seeds on each of the simulation chunks (useful for reproducing parallel simulations). Specifically, for k in 1:chunks, seeds are set as set.seed(seeds[k], kind = "Mersenne-Twister") in each chunk. Defaults to NULL (no seed setting is done).

...

optional arguments to be passed to f or to foreach (for example, .export to export global variables or other functions to the foreach environment).

Details

It is possible to have a progress bar if int_sph_MC is wrapped with progressr::with_progress or if progressr::handlers(global = TRUE) is invoked (once) by the user. See the examples below. The progress bar is updated with the number of finished chunks.

Examples

## Sequential simulation

# Vectorized functions to be integrated
x1 <- function(x) x[, 1]
quad <- function(x, a = 0) a + rowSums(x^4)

# Approximate \int_{S^{p-1}} x_1 dx = 0
int_sph_MC(f = x1, p = 3, M = 1e4, chunks = 2)

# Approximate \int_{S^{p-1}} (a + \sum_i x_i^4) dx
int_sph_MC(f = quad, p = 2, M = 1e4, a = 0, chunks = 2)

# Compare with Gauss--Legendre integration on S^2
th_k <- Gauss_Legen_nodes(a = 0, b = 2 * pi, N = 40)
w_k <- Gauss_Legen_weights(a = 0, b = 2 * pi, N = 40)
sum(w_k * quad(cbind(cos(th_k), sin(th_k)), a = 1))

## Parallel simulation with a progress bar
# \donttest{
# Define a progress bar
require(progress)
require(progressr)
handlers(handler_progress(
  format = paste("(:spin) [:bar] :percent Iter: :current/:total Rate:",
                 ":tick_rate iter/sec ETA: :eta Elapsed: :elapsedfull"),
  clear = FALSE))
# Call int_sph_MC() within with_progress()
with_progress(int_sph_MC(f = x1, p = 3, cores = 2, M = 1e5, chunks = 100))

# Instead of using with_progress() each time, it is more practical to run
# handlers(global = TRUE)
# once to activate progress bars in your R session
# }