%matplotlib notebook
import math
import numpy as np
import matplotlib.pyplot as plt

Sampling and integration in 1D

Here is a function I might like to integrate.

def f(x):
    return 1 - np.sqrt(1 - x**4)

plt.figure()
xs = np.linspace(0, 1, 501)
plt.plot(xs, f(xs))

[<matplotlib.lines.Line2D at 0x123e3fc18>]

The integral of this function doesn't have a closed form. Wolfram Alpha can integrate it in terms of the Elliptic Integral and gives the approximate value 0.125981 for the integral from 0 to 1.

The simplest Monte Carlo approach is to use uniform samples over the desired integration interval.

def estimate_simple(f, N):
    """Integrate f from 0 to 1 using N uniform random samples."""
    xs = np.random.random_sample(N)
    return np.mean(f(xs))

estimate_simple(f, 1000000) - 0.125981

-0.00025539867922452775

I'm going to package this up in a "sampling strategy" that has the function, the code to generate samples, and the pdf of the samples, organized into a standard 3-method interface.

strategy1 = {
    'eval': lambda x: 1 - np.sqrt(1 - x**4),
    'sample': lambda xi: xi,
    'pdf': lambda x: np.ones(x.shape)
}

def estimate(strategy, N):
    """Integrate a function using N independent samples drawn from a provided pdf."""
    xis = np.random.random_sample(N)
    xs = strategy['sample'](xis)
    return np.mean(strategy['eval'](xs) / strategy['pdf'](xs))

estimate(strategy1, 100) - 0.125981

-0.024124005829066306

def stdev_of(fn, M):
    """Compute the standard deviation of repeated calls to a function."""
    results = [fn() for i in range(M)]
    return np.std(results)

stdev_of(lambda: estimate(strategy1, 10000), 100)

0.00019269334701036699

strategy2 = {
    'eval': lambda x: 1 - np.sqrt(1 - x**4),
    'sample': lambda xi: xi ** (1/5),
    'pdf': lambda x: 5 * x**4
}

def histcheck(strategy, N):
    xis = np.random.random_sample(N)
    xs = strategy['sample'](xis)
    ps, edges, patches = plt.hist(xs, math.ceil(math.sqrt(N)/10), density=True)
    plt.plot(edges, strategy['pdf'](edges), 'k')

plt.figure()
histcheck(strategy2, 1000000)

estimate(strategy2, 100) - 0.125981

0.0010083955354215912

stdev_of(lambda: estimate(strategy2, 10000), 100)

0.00018072102690918401

Sampling uniformly from unit disc

To generate samples in the unit disc it is convenient to generate them in polar coordinate form. Unfortunately uniformly generating r and theta does not do it:

rs = np.random.random_sample(10000)
thetas = 2 * np.pi * np.random.random_sample(10000)
plt.figure()
plt.plot(rs * np.cos(thetas), rs * np.sin(thetas), '.', ms=1)
plt.axis('equal')

rs = np.sqrt(np.random.random_sample(10000))
plt.figure()
plt.plot(rs * np.cos(thetas), rs * np.sin(thetas), '.', ms=1)
plt.axis('equal')