Monte Carlo integration is a technique for numerical integration using random numbers. It can be useful, for example, when evaluating integrals over domains of high dimension. Meshing in high dimensions suffers from the [curse of dimensionality](https://en.wikipedia.org/wiki/Curse_of_dimensionality): 100 evenly spaced sample points suffice to sample the unit interval $[0,1]$ with no more than 0.01 distance between points, whereas sampling of the 10-dimensional unit hypercube $[0,1]^{10}$ with a lattice that has a spacing of 0.01 between adjacent points would require $10^{20}$ sample points.

Let $\Omega \subset \mathbb R^n$ and let $f : \Omega \to \mathbb R$ be piecewise continuous. 
Let $p$ be such a probability density function on $\Omega$ that $p(x) > 0$ for all $x \in \Omega$.
The Monte Carlo approach to approximate the integral $\int_\Omega f(x) dx$
works as follows: 

1. Generate a large number of independent samples $x_1,\dots,x_N \in \Omega$
from the probablity distribution with the density $p$
2. Compute the quantity 

$$
I_N = \frac 1 N \sum_{i=1}^N \frac{f(x_i)}{p(x_i)}.
$$

The [law of large numbers](https://en.wikipedia.org/wiki/Law_of_large_numbers) implies that 

$$
I_N \to \int_\Omega f(x) dx \quad \text{as $n \to \infty$}.
$$

Choosing $p$ cleverly is the basic idea behind [importance sampling](https://en.wikipedia.org/wiki/Importance_sampling).

[Buffon's needle](https://en.wikipedia.org/wiki/Buffon%27s_needle_problem) experiment is a classical pedagogical example of Monte Carlo integration. 
Suppose we have a floor made of parallel strips of wood, each the same width $t$, and we drop a needle of length $l < t$ onto the floor. Buffon showed that the  probability $P$ that the needle will lie across a line between two strips is

\begin{equation}\tag{1}
P=\frac{2}{\pi}\frac{l}{t}.
\end{equation}

Let $s$ be the distance from the center of the needle to the closest parallel line, and let $\theta$ be the acute angle between the needle and one of the parallel lines. 
Here $s$ and $\theta$ are random variables with uniform distributions over $[0, t/2]$ and $[0,\pi/2]$, respectively. We write $p(s,\theta)$ for their joint probability density function and $A$ for the event that the needle lies across a line between two strips. Then 

$$
P = \int_A p(s,\theta) ds d\theta.
$$

Writing $\Omega = [0, t/2] \times [0,\pi/2]$ and $f(x) = 1_A(x) p(x)$
where $x = (s,\theta)$ and

$$
1_A(x) = 
\begin{cases}
1 & x \in A,
\\
0 & x \notin A,
\end{cases}
$$

we have, using the notation above, that $I_N \to P$ as $N \to \infty$.
Observe that, in this case,

$$
I_N = \frac 1 N \sum_{i=1}^N 1_A(x_i),
$$

and evaluating $I_N$ boils down to counting the needles that lie across a line between two strips.

Take $l = 5/6$ and $t = 1$. The goal of this homework is to approximate 
$\pi$ via the formula 

$$
\pi=\frac{2l}{t P} \approx \frac{2l}{t I_N},
$$

that follows from (1). 

Let us give some remarks on the history of Monte Carlo integration. Formula (1) was first derived in 

> Buffon, Georges L. L., comte de. _Histoire naturelle, générale et particulière, Supplément 4_. Imprimerie royale, Paris, 1777. (scan in [Google Books](https://books.google.fi/books?id=AjhYD1vsVAIC&hl=fi&pg=PA100#v=onepage&q&f=false))

Now my French is not very strong, but as claimed [here](https://en.wikipedia.org/wiki/Buffon's_needle_problem), it seems that approximating $\pi$ was not the original motivation for Buffon's question. You can have a look at his book, see pp. 100-104, and while you are at it, you can also try figure out if there is an error in his derivation, as claimed [here](https://mathworld.wolfram.com/Buffon-LaplaceNeedleProblem.html). 

The idea of using Buffon's formula to design a method for approximating the number $\pi$ goes back at least to Laplace. _"Si l'on projette un grand nombre de fois ce cylindre [...] ce qui fera connaître la valeur de la circonférence $2 \pi$"_, see p. 360 of

> Laplace, Pierre S., marquis de. _Théorie analytique des probabilités_. Veuve Courcier, Paris, 1812. (scan in [Internet Archive](https://archive.org/details/thorieanalytiqu01laplgoog/page/n464/mode/2up))

The first computerized Monte Carlo simulations were run on [ENIAC](https://en.wikipedia.org/wiki/ENIAC) in 1948 by a team including John and Klara von Neumann and Nick Metropolis. It can be argued that the simulations were also the first code written in the modern paradigm, associated with the "stored program concept," ever to be executed, see 

> Haigh, Thomas, Priestley, Mark, and Rope, Crispin. _Los Alamos Bets on ENIAC: Nuclear Monte Carlo Simulations, 1947-1948_. IEEE Annals of the History of Computing 36, no. 3, 42-63, 2014. <https://doi.org/10.1109/MAHC.2014.40> (in [Helka](https://helka.helsinki.fi/permalink/358UOH_INST/qn0n39/cdi_ieee_primary_6880250))



```python
import numpy as np
rng = np.random.default_rng()

t = 1
l = 5/6

def sample():
    '''Returns s and theta generated using the random number generator rng'''
    # Draw samples from uniform distributions using the function rng.uniform, see
    # https://numpy.org/doc/stable/reference/random/generated/numpy.random.Generator.uniform.html
    raise NotImplementedError() # a placeholder, your implementation goes here
    
def intersects(s, theta):
    '''Returns True iff the needle lies across a line between two strips'''
    raise NotImplementedError() # a placeholder, your implementation goes here

# We will plot I_N for several N so let's save every nth approximation
def I(n, K):
    '''Return I_n, I_{2n}, ..., I_{Kn}'''
    out = np.zeros(K)
    raise NotImplementedError() # a placeholder, your implementation goes here
    return out
```


```python
import matplotlib.pyplot as plt
```


```python
# Plot one sample
s, theta = sample()
c = np.array([0, s]) # Center of needle
d = np.array([np.cos(theta), np.sin(theta)]) # Direction of needle
# End points of needle
end1 = c - l/2*d
end2 = c + l/2*d
ends = np.stack((end1, end2))
xs = ends[:,0]
ys = ends[:,1]

plt.plot(xs, ys, 'r') # needle in red
plt.plot([-1,1],[0, 0], 'b') # closest line in blue
ax = plt.gca()
ax.set_ylim(-0.5, 1)
ax.set_aspect(1)
print(f'Needle intersects the closest line: {intersects(s, theta)}')
```


```python
# Plot convergence to pi
n = 100
K = 40
Is = I(n, K)
Ns = n*np.arange(1,K+1)
plt.plot(Ns, 2*l/(t*Is), 'b') # approximation in blue
plt.plot([n,n*K],[np.pi, np.pi],'r') # pi in red
ax = plt.gca()
ax.set_ylim(2.4, 3.6);
```

**How to hand in your solution**

1. Run the whole notebook by choosing _Restart Kernel and Run All Cells_ in the _Run_ menu
    - Alternatively you can click the ⏩️ icon in the toolbar
2. Click the link below to check that the piece of code containing your solution was uploaded to pastebin
    - If you have changed the order of cells in the notebook, you may need to change the number in the below cell to the one in the left margin of the cell containing your solution
3. Copy the link and submit it in Moodle
    - You can copy the link easily by right-clicking it and choosing _Copy Output to Clipboard_


```python
# Upload the code in the first input cell to pastebin
%pastebin 1
```

0da74dd9bb865551155a71f1e1d9af972f5fc7ab
ipynb
Jupyter Notebook
integration/homework.ipynb
uh-comp-methods1/notebooks
8c8586094f7fba06053f9dd5ba370ebd87389645
CC-BY-4.0
2022-01-11T13:08:56.000Z
2022-01-11T13:09:47.000Z
integration/homework.ipynb
uh-comp-methods1/notebooks
8c8586094f7fba06053f9dd5ba370ebd87389645
CC-BY-4.0
integration/homework.ipynb
uh-comp-methods1/notebooks
8c8586094f7fba06053f9dd5ba370ebd87389645
CC-BY-4.0
2022-01-18T09:48:10.000Z
2022-01-24T18:03:41.000Z
Qwen/Qwen-72B
1. YES
2. YES