Monte-Carlo Integration

For uniformly sampled points $x_i$ in a domain of volume $\mathrm{\Omega }$, a basic Monte-Carlo estimator is

$${\hat{I}}_N=\frac{\mathrm{\Omega }}{N}\sum _{i=1}^{N}f(x_i).$$

This example estimates ${\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}{e}^{-x^2}dx$ by sampling the finite interval [-10, 10]. @accelerate frees the non-returned sample arrays after their final use. CUDA may also fuse eligible broadcasts.

# examples/integrate.jl
using cuNumeric

integrand(x) = exp(-x^2)

@accelerate function monte_carlo(N, x_max)
    Ω = 2 * x_max
    raw_samples = cuNumeric.rand(N)
    samples = @. Ω * raw_samples - x_max # first kernel
    return (Ω / N) * sum(integrand.(samples)) # second (integrand) + third kernels (sum)
end

const N = 1_000_000
const x_max = 10.0f0
estimate = monte_carlo(N, x_max)
println("Monte-Carlo estimate: $(fetch(estimate))")
println("Analytical value: $(sqrt(pi))")

The result is a CNScalar backed by a 0D NDArray, which keeps the reduction asynchronous. Use fetch(estimate) only when a Julia scalar is required.

This can also be done with mapreduce which fuses everything into a single kernel.

const Ω = 2 * x_max
mc = (x) -> integrand(Ω*x - x_max)
estimate = (Ω/N) * mapreduce(mc, +, cuNumeric.rand(N); init = 0.0f0) # only kernel