---
author:
- "by [Eugeniy E. Mikhailov](http://physics.wm.edu/~evmik/) and [Greg
  Bentsen](https://physics.wm.edu/~gbentsen/)"
title: Monte-Carlo numerical integration
---

## Logistics and Agenda

-   Monday: Monte-Carlo integration
-   Tuesday: root finding and integration gotchas
-   Friday: mock in-class exam

## Midterm 1:

Do it at home part

-   Midterm write up and code is due Sunday October 4th 11:59pm
    -   **there is no late submission option**
-   You may **not** use AI to write the report.
    -   Report should be done by you and only one per group is needed

In class part

-   on Friday October 2nd: mock exam (no grades and you do not even have
    to do it)
    -   check if your computers are battle worthy for off-line use
    -   see a typical problem set
-   on Monday October 5th, with computers but no AI and no Internet
    connection
    -   we will give you modified assignment which would be similar
    -   it will require modification or reuse of prepared code

# Monte Carlo Integration

## Toy example - area of the pond

::: columns
::: {.column width="55%"}
![Pond and random hits illustration](./pics/pond_area.svg){width="100%"}
:::

::: {.column width="45% "}
::: definition
$$
A_{pond}=\frac{N_{inside}}{N_{total}} A_{box}
$$ where $$
A_{box} = (b_x-a_x)  (b_y-a_y)
$$
:::

-   Points must be `\alert{uniformly}`{=tex} and randomly distributed
    across the area.
-   The smaller the enclosing box the better it is.
:::
:::

## Naive Monte Carlo integration

::: columns
::: {.column width="55%"}
![Naive Monte-Carlo setup](./pics/monte_carlo_integr.svg){width="100%"}
:::

::: {.column width="45% "}
::: definition
$$
\int_{a_x}^{b_x} f(x) dx =\frac{N_{inside}}{N_{total}} A_{box}
$$ where $$
A_{box} = (b_x-a_x)  (b_y)
$$
:::

-   Points must be `\alert{uniformly}`{=tex} and randomly distributed
    across the area.
-   The smaller the enclosing box the better it is. So
    $max(f(x)) \to b_y$
:::
:::

## Monte Carlo integration derived

::: columns
::: {.column width="55%"}
![Improved Monte-Carlo
setup](./pics/monte_carlo_integr_demistified.svg){width="100%"}
:::

::: {.column width="45% "}
Notice that if we choose a small stripe around the bin value $x_b$, then
subset of points in that stripe gives an estimate for $f(x_b)$.

Thus why bother spreading points around area?

::: definition
Let's chose a uniform random distribution of points ${x_i}$ inside
$[a_x, b_x]$

$$
\int_{a_x}^{b_x} f(x) dx \approx \frac{b_x - a_x}{N}\sum_{i=1}^N f(x_i)
$$
:::
:::
:::

# Numerical Integration methods error estimates

## Error estimate for Monte-Carlo method

::: definition
**Monte Carlo method algorithmic error $(E)$**

$$
E={\mathbf{O}}\left( (b_x -a_x) {\frac{\sigma_f} {\sqrt{N}} } \right)
$$

where $\sigma_f$ is the estimate of the function standard deviation
given by

```{=tex}
\begin{aligned}
            \sigma_f &=& \sqrt{<f^2>-<f>^2} \\
            <f>   &=& \lim_{M \to \inf} \frac{1}{M}\sum_{i=1}^M f(x_i) \\
            <f^2> &=& \lim_{M \to \inf} \frac{1}{M}\sum_{i=1}^M f^2(x_i)
\end{aligned}
```
:::

On the first glance it is horendous method in comaparison even to the
box method where $E \sim 1/N$

## Monte-Carlo 1D example

::: columns
::: {.column width="45% "}
``` python
import numpy as np
def monte_carlo(func, a, b, N=1_000_000):
    """Monte Carlo Integration"""
    xi = np.random.uniform(a, b, N)
    f_xi = func(xi)
    return np.mean(f_xi) * (b - a)

if __name__ == "__main__":
    import matplotlib.pyplot as plt
    def f(x):
        return x * x
    a = 0
    b = 1
    Np = 100  # hom many points to plot
    N = np.logspace(1, 6, Np, dtype=np.int64)  # used in Monte-Carlo
    integral_for_N = np.zeros(Np)
    for i, n in enumerate(N):
        integral_for_N[i] = monte_carlo(f, a, b, n)

    plt.loglog(N, np.abs(integral_for_N - 1 / 3))
    plt.grid(True)
    plt.xlabel("N")
    plt.ylabel("Deviation from true value")
    plt.savefig("monte_carlo_deviation.png")
```
:::

::: {.column width="55%"}
![Scaling of Monte-Carlo algorithmic
error](./code/monte_carlo_deviation.png){width="100%"}
:::
:::

# Multidimensional integration

## Multidimensional integration with interval splitting

$$
        \int_{a_x}^{b_x} 
        \int_{a_y}^{b_y} 
        f(x,y) ~ dx ~ dy  =
        \int_{a_x}^{b_x} dx
        \int_{a_y}^{b_y} dy
        ~ f(x,y) 
$$

Note that last integral is the function of only $x$

```{=tex}
\begin{aligned}
        \int_{a_y}^{b_y} dy ~ f(x,y)  = F(x)
    
\end{aligned}
```
```{=tex}
\begin{aligned}
        \int_{a_x}^{b_x} 
        \int_{a_y}^{b_y} 
        f(x,y) ~ dx ~ dy  =
        \int_{a_x}^{b_x} dx
        ~ F(x) 
    
\end{aligned}
```
Thus we replaced multidimensional integral as consequent series of
single dimension integrals, which we already know how to do.

3D case would look like this

$$
        \int_{a_x}^{b_x} 
        \int_{a_y}^{b_y} 
        \int_{a_z}^{b_z} 
        f(x,y,z) ~ dx ~ dy ~ dz =
        \int_{a_x}^{b_x} dx
        \int_{a_y}^{b_y} dy
        \int_{a_z}^{b_z} dz 
        ~ f(x,y,z) 
$$

## Multidimensional integration with Monte Carlo

Note that if we would like to split integration region by $N$ points in
every of $D$ dimensions, then evaluation time grows $\sim N^D$, which
renders Rectangle, Trapezoidal, Simpson, and alike methods useless.

Monte Carlo method is a notable exception, it looks very simple even for
multidimensional case and maintains the same $\sim N$ evaluation time.

3D case would look like this

$$
        \int_{a_x}^{b_x} dx
        \int_{a_y}^{b_y} dy
        \int_{a_z}^{b_z} dz 
        ~ f(x,y,z) 
        \approx
        \frac{ 
            (b_x - a_x)
            (b_y - a_y)
            (b_z - a_z)
        }{N}\sum_{i=1}^N f(x_i, y_i, z_i)
$$

**A general case**

$$
        \int_{V} dV
        ~ f(\vec{x}) 
        =
        \int_{V} dx_1 dx_2 dx_3 ... dx_D
        ~ f(\vec{x}) 
        \approx
        \frac{V}{N}
        \sum_{i=1}^N f(\vec{x}_i)
$$

where V is the multidimensional volume , ${\vec{x}_i}$ randomly and
uniformly distributed points in the volume V.

# Python library Monte-Carlo method

Quasi-Monte Carlo quadrature `qmc_quad` from `scipy.integrate` package

``` python
qmc_quad(func, a, b, *, n_estimates=8, n_points=1024, qrng=None, log=False)
```

it has many more parameters then our simple function

``` python
from scipy.integrate import qmc_quad
def f(x): return x*x
res = qmc_quad(f, 0, 1, n_points=10000)

In [28]: res
Out[28]: QMCQuadResult(integral=np.float64(0.3333562442002344), standard_error=np.float64(2.1699395303402513e-05))

In [29]: res.integral
Out[29]: np.float64(0.3333562442002344
```

Here is example of 2D integral, which finds the area of the disk with
radius 1, i.e. the final result should be $\pi$

``` python
In [104]: from scipy.integrate import qmc_quad

In [105]: def disk(rv): return np.int64(np.sum(rv*rv, axis=0) <=1)

In [106]: res=qmc_quad(disk, [-1, -1], [1,1], n_points=10000)

In [107]: res
Out[107]: QMCQuadResult(integral=np.float64(3.141000000000001), standard_error=np.float64(0.000994269293788858))

In [108]: res.integral
Out[108]: np.float64(3.141000000000001)
```