In [1]:

import numpy as np

Simple Monte Carlo estimators¶

This notebook shows some very simple examples of Monte Carlo integration we played with in lecture. If you're looking at the HTML export of this notebook and would like to play with it interactively, the source file is here.

Here is one of the random number generators in the Python library.

In [2]:

np.random.random_sample()

Out[2]:

0.09396070581276861

The expected value of the RNG output is 0.5. This is easier to see if we average multiple trials:

In [3]:

np.mean(np.random.random_sample(10000000))

Out[3]:

np.float64(0.4999689107001966)

Let's say we want to integrate the function f(x):

In [4]:

def f(x):
    return np.sin(x)

Integral of $f(x)$ from 0 to 1 is $1 - \cos(1)$.

In [5]:

1 - np.cos(1)

Out[5]:

np.float64(0.45969769413186023)

In [6]:

np.mean(f(np.random.random_sample(10000000)))

Out[6]:

np.float64(0.4595575960516077)

So $f(x)$ is an estimator for the integral of $f(x)$ over the unit interval, using unit random numbers. What if we want to integrate from $0$ to $\pi/2$, where we know the integral evaluates to $1.0$? We need samples on the interval $[0,\pi/2]$, which we can get by scaling up our unit random numbers.

In [7]:

np.mean(f(np.random.random_sample(1000000) * np.pi/2))

Out[7]:

np.float64(0.6366494086959851)

Hmmm, this is not the right answer. What should the estimator be?