int_sph_MC function

Monte Carlo integration of functions on the (hyper)sphere

Monte Carlo approximation of the integral [REMOVE_ME]$\int_{S^{p-1}}f(x)\,\mathrm{d}x\int_{S^{p-1}} f(x) dx [REMOVE_ME_2]$

of a function $f:S^{p-1} \rightarrow R$ defined on the (hyper)sphere c("$\n$", "$S^{p-1}:={x\\in R^p:||x||=1}$"), $p\ge 2$.

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

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).

Returns

A scalar with the approximate integral.

Description

Monte Carlo approximation of the integral

of a function $f:S^{p-1} \rightarrow R$ defined on the (hyper)sphere c("$\n$", "$S^{p-1}:={x\\in R^p:||x||=1}$"), $p\ge 2$.

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

# 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