```
%matplotlib notebook
import math
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
```

# 1 Some Monte Carlo integration examples

This notebook shows some examples of Monte Carlo integration which we’ll look at in class. Some are 1D integrals on the unit integral; others are 2D integrals on the unit hemisphere, which lead directly to algorithms we use for illumination.

If you are looking at the HTML export of this notebook, the source file is here: monte-carlo.ipynb and you can download it and run it in Jupyter Notebook if you have it. (If you’d like to get it, installing Anaconda is a quick way.)

## 1.1 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()= np.linspace(0, 1, 501)
xs ; plt.plot(xs, f(xs))
```

`<IPython.core.display.Javascript object>`

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."""
= np.random.random_sample(N)
xs return np.mean(f(xs))
```

I can try this estimator with different numbers of samples to see how the accuracy improves:

`10000000) - 0.125981 estimate_simple(f, `

`-4.7454293878868015e-06`

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."""
= np.random.random_sample(N)
xis = strategy['sample'](xis)
xs return np.mean(strategy['eval'](xs) / strategy['pdf'](xs))
```

I can do the same experiment as before using this interface.

`1000000) - 0.125981 estimate(strategy1, `

`-1.717783785970539e-05`

Let’s more scientifically compute the standard deviation of this estimator as a function of the number of samples.

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

`lambda: estimate(strategy1, 1000000), 100) stdev_of(`

`0.0001795481760075688`

```
= np.int32(np.logspace(2, 5, 7))
Ns = [stdev_of(lambda: estimate(strategy1, N), 100) for N in Ns]
stdevs
plt.figure()'-')
plt.loglog(Ns, stdevs, '# samples')
plt.xlabel('standard deviation');
plt.ylabel('equal'); plt.axis(
```

`<IPython.core.display.Javascript object>`

That shows a very clear trend: standard deviation is proportional to the square root of the number of samples.

Here is another sampling stragegy, using the probability density .

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

The `sample`

method is using the inverse-cdf
algorithm. The desired pdf is ;
the corresponding cumulative distribution function (cdf) is

and the inverse cdf is

First let’s verify that this sampling method really generates
the pdf that it claims to. The following function computes many
samples drawn from `sample`

and plots a histogram
showing the density of those samples against the value returned by
the `pdf`

method:

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

```
plt.figure()100000) histcheck(strategy2,
```

`<IPython.core.display.Javascript object>`

Good, that shows that the `pdf`

method really does
report the probability with which the `sample`

method
generates a sample near a given value.

`1000000) - 0.125981 estimate(strategy2, `

`-2.196604269441571e-05`

`lambda: estimate(strategy2, 1000000), 100) stdev_of(`

`1.993623186073966e-05`

Now let’s compare the standard deviation of this estimator to the other one, for various numbers of samples.

```
= np.int32(np.logspace(2, 5, 7))
Ns = [stdev_of(lambda: estimate(strategy1, N), 100) for N in Ns]
stdevs1 = [stdev_of(lambda: estimate(strategy2, N), 100) for N in Ns]
stdevs2
plt.figure()
plt.loglog(Ns, stdevs1, Ns, stdevs2)'# samples')
plt.xlabel('standard deviation');
plt.ylabel('equal'); plt.axis(
```

`<IPython.core.display.Javascript object>`

This shows that second strategy outperforms the first one by a constant factor of about 10.

## 1.2 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 and does not do it:

```
def plot_circle(r, c='k'):
* np.cos(np.linspace(0,2*np.pi,121)), r * np.sin(np.linspace(0,2*np.pi,121)), c) plt.plot(r
```

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

`<IPython.core.display.Javascript object>`

These samples are clumped towards the middle because the radii are uniformly distributed, so there are roughly equal numbers of points with and . But the areas on the plane corresponding to these two intervals in are very different.

```
0.1)
plot_circle(0.9) plot_circle(
```

The area occupied by the interval between and is proportional to , so to generate points uniformly we want . The integral of this is , so the normalized pdf is with the cdf . This means if we have a random number that’s uniform on the unit interval, we can compute and then .

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

`<IPython.core.display.Javascript object>`