# Physics 256
## Rejection Sampling





## Last Time

- Generation and testing of pseudorandom numbers
- Tower sampling

## Today

- Exponential distribution
- Box-Muller for Gaussian random numbers
- Rejection Sampling


## Setting up the Notebook


```python
import matplotlib.pyplot as plt
import numpy as np
%matplotlib inline
plt.style.use('notebook');
%config InlineBackend.figure_format = 'retina'
colors = ["#2078B5", "#FF7F0F", "#2CA12C", "#D72827", "#9467BE", "#8C574B",
            "#E478C2", "#808080", "#BCBE20", "#17BED0", "#AEC8E9", "#FFBC79", 
            "#98E08B", "#FF9896", "#C6B1D6", "#C59D94", "#F8B7D3", "#C8C8C8", 
           "#DCDC8E", "#9EDAE6"]
```

## Sampling Continous Distributions

Last time we saw that we can sample from any continuous probability distribution $p(y)$ using only uniformly distributed random numbers $x$:

\begin{equation}
P(y) = \int_{-\infty}^y p(y') dy' = \int_{-\infty}^{y} \frac{dx}{dy'}  p_u(x) dy' = \int_0^y \frac{dx}{dy'} dy' = x(y)
\end{equation}

So $y = P^{-1}(x)$ for $x \in \mathcal{U}_{[0,1]}$.


### Example 2
Consider sampling from the exponential distribution for $y\ge 0$:

\begin{equation}
p(y) = \lambda \mathrm{e}^{-\lambda y}.
\end{equation}

As before:

\begin{equation}
x(y) =\int_0^y p(y') dy' = \lambda \int_0^y \mathrm{e}^{-\lambda y'}dy' = \left.-\mathrm{e}^{-\lambda y'} \right \rvert_0^y = 1 - \mathrm{e}^{-\lambda y} .
\end{equation}

Solving for $x$:

\begin{align*}
\mathrm{e}^{-\lambda y} &= 1-x \newline
-\lambda y &= \ln (1-x) \newline
y &= -\frac{1}{\lambda} \ln (1-x) 
\end{align*}

which gives $y \in [0,\infty)$ for $x \in [0,1)$.

<br />
<div class="span alert alert-danger">
Note: we need to be careful with our uniform random number $x$ as it can formally take on the end point value $x=1$. 
</div>

<div class="span alert alert-success">
<h2>Programming challenge </h2>
Generate $N=1000$ exponentially distributed random numbers with $\lambda = 0.5$ and check their distribution using a histogram.
</div>


```python
def p(y,λ):
    return λ*np.exp(-λ*y)

N = 100000
λ = 0.5
y = np.linspace(0,100,N)
plt.plot(y,p(y,λ),color=colors[0], label=r'$%3.1f\mathrm{e}^{-%3.1fy}$'%(λ,λ))

# sample y from a uniform x
x = np.random.random(N)
sampled_y = -(1/λ)*np.log(1.0-x)

plt.hist(sampled_y, bins=100, normed=True, ec='w', label='sampled', fc=colors[0], alpha=0.2);
plt.xlim(0,10);
plt.xlabel('y')
plt.ylabel('p(y)')
plt.title('Exponential Distribution')
plt.legend(loc='upper right')
```

## Gaussian Distributed Random Numbers

Consider the normal or Gaussian probability distribution with mean $\mu=0$:

\begin{equation}
p_n(y) = \frac{1}{\sqrt{2\pi\sigma^2}} \mathrm{e}^{-y^2/2\sigma^2}.
\end{equation}

We have already seen the cumulative distribution function: 

\begin{align*}
P_n(y) &= \int_{-\infty}^{y} p_n(y') \newline
&=\frac{1}{2}\left[1 + \mathrm{erf}\left(\frac{y}{\sqrt{2} \sigma}\right) \right]
\end{align*}

which cannot be analytically inverted.  However, we can use the same trick that is used to evaluate Gaussian integrals, i.e. extending our calculation to higher dimensions.  In this case, that means considering the multivariate probability distribution:

\begin{equation}
p_n(x,y) = p_n(x)p_n(y) = \frac{1}{2\pi\sigma^2} \mathrm{e}^{-(x^2+y^2)/2\sigma^2}.
\end{equation}

Let us transform to polar coordinates defined by:

\begin{align*}
r &= \sqrt{x^2+y^2} \newline
\theta &= \tan^{-1} \left(\frac{y}{x}\right)
\end{align*}

thus:

\begin{equation}
p_n(x,y)dx dy = \frac{1}{2\pi\sigma^2} \mathrm{e}^{-r^2/2\sigma^2} r dr d\theta .
\end{equation}

Defining

\begin{equation}
\rho = \frac{r^2}{2\sigma^2}  \Rightarrow d\rho = \frac{r dr}{\sigma^2}
\end{equation}

we have:

\begin{align*}
x &= \sqrt{2\rho}\sigma \cos\theta \newline
y &= \sqrt{2\rho}\sigma \sin\theta \newline
\end{align*}

and

\begin{equation}
p_n(x,y)dx dy = \frac{1}{2\pi} \mathrm{e}^{-\rho} d\rho d\theta .
\end{equation}

We know how to sample both scaled uniform random numbers $\theta \in \mathcal{U}_{[0,2\pi]}$ and exponentially distributed random numbers $\rho \in \mathrm{e}^{-\rho}$.  Therefore we can write:

\begin{align*}
x &= \sqrt{-2\ln(1-x_1)}\sigma \cos(2\pi x_2) \newline
y &= \sqrt{-2\ln(1-x_1)}\sigma \sin(2\pi x_2)
\end{align*}

where $x_1,x_2 \in \mathcal{U}_{[0,1)}$, i.e. 2 uniformly distributed random numbers can be used to get two Gaussian distributed random numbers.  This is known as the **Box-Muller** method.


```python
from scipy.constants import pi as π
N = 100000

# our uniform random numbers
x1 = np.random.random(N)
x2 = np.random.random(N)

# generate the Box-Muller values
x = np.sqrt(-2.0*np.log(1-x1))*np.cos(2.0*π*x2)
y = np.sqrt(-2.0*np.log(1-x1))*np.sin(2.0*π*x2)

# combine them into 1 array
r = np.hstack([x,y])

# produce a plot comparing with the actual distribution
px = np.linspace(-4,4,1000)
plt.plot(px,np.exp(-px**2/2)/np.sqrt(2*π), color=colors[0], label=r'$\frac{1}{\sqrt{2\pi}}\mathrm{e}^{-y^2/2}$')
plt.hist(r, bins=100, normed=True, ec='w', label='sampled',fc=colors[0], alpha=0.2)
plt.xlabel('y')
plt.ylabel('p(y)')
plt.xlim(-4,4)
plt.title('Box-Muller')
plt.legend(loc='upper left')
```

## Rejection Sampling

What do we do if we can't invert the cumulative distribution function and we are out of tricks?  The technique, known as **von Neumann rejection** is similar to how we computed $\pi$ by throwing stones into a pond.  

Consider some probability distribution function $p(x)$ whose cumulative distribution function $P(x)$ cannot be inverted.  Suppose we want to sample this function over the finite range $[x_{\rm min},x_{\rm max}]$ where $p_{\rm max}$ is the maximal value of $p(x)$ on this interval.  Then we can:

1. Generate a sequence of $N$ uniformly distributed random numbers $x_i \in \mathcal{U}_{[x_{\rm min},x_{\rm max}]}$
2. Generate a second sequence of $N$ uniformly distribted random numbers $y_i \in \mathcal{U}_{[0,p_{\rm max}]}$
3. Keep only those elements in the first sequence $x_i$ that have $y_i < p(x_i)$

The resulting set of $x_j$ values will be distributed according to $p(x)$.   This works by throwing away unlikely results.

<br />
<div class="span alert alert-warning">
Note: Our probability distribution must be bounded and have a finite range. 
</div>

Easiest to understand via an example.

### Example
Sample $N=10^4$ random numbers from the function:

\begin{equation}
f(x) = \frac{1}{\sqrt{\cosh(x)}}
\end{equation}

on the region $[-5,5]$.

#### Step 1: at present, $f(x)$ it not normalized so we need to turn it into a probability distribution.


```python
from scipy import integrate
xmin,xmax = -5,5

f = lambda x: 1.0/np.sqrt(np.cosh(x))
A = 1.0/integrate.quad(f, xmin, xmax)[0]
print(A)

p = lambda x: A*f(x)
```

    0.2092148976418844


#### Step 2: determine the maximual value of $p(x)$ on the interval


```python
# Let's start by plotting
px = np.linspace(xmin,xmax,10000)
plt.plot(px,p(px), color=colors[0], label=r'$\frac{%4.2f}{\sqrt{\cosh(x)}}$'%A)
plt.xlabel('x')
plt.ylabel('p(x)')
plt.xlim(-5,5);
plt.legend()
```

Thus, the maxima will appear at $x=0$.  We could also use [``scipy.optimize.minimize``](https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.minimize.html#scipy.optimize.minimize) to minimize $-p(x)$


```python
print('p(0) = ',p(0))

from scipy.optimize import minimize
pmax = -minimize(lambda x: -p(x),-5).fun
print('max p(x) = ',pmax)
```

    p(0) =  0.209214897642
    max p(x) =  0.20921489764188417


#### Step 3: Generate the two lists of random numbers


```python
N = 10**6
x = xmin + (xmax-xmin)*np.random.random(N)
y = p(0)*np.random.random(N)
```

#### Step 4: perform the rejection


```python
#accepted = np.array([])
#for i in range(N):
#    if y[i] < p(x[i]):
#        accepted = np.append(accepted,x[i])
        
accepted = x[y < p(x)]
```

#### Check that it works!


```python
plt.plot(px,p(px),color=colors[0], label=r'$\frac{%4.2f}{\sqrt{\cosh(x)}}$'%A)
plt.hist(accepted, bins=100, normed=True, ec='w', label='sampled', fc=colors[0], alpha=0.2)
plt.xlabel('x')
plt.ylabel('p(x)')
plt.xlim(-5,5);
plt.legend()
```

#### Let's see what the process actually looks like


```python
x_accepted = x[y < p(x)]
y_accepted = y[y < p(x)]

x_rejected = x[y >= p(x)]
y_rejected= y[y >= p(x)]

plt.axhline(y=pmax, color='gray', lw=2, label='p(0)')
plt.plot(px,p(px),color=colors[0], label=r'$\frac{%4.2f}{\sqrt{\cosh(x)}}$'%A)

plt.plot(x_accepted,y_accepted,'o', mfc='None', mec=colors[0], ms=3, label='accepted')
plt.plot(x_rejected,y_rejected,'x', mfc='None', mec=colors[1], ms=3, label='rejected')
plt.xlabel('x')
plt.ylabel('p(x)')
plt.xlim(-5,5);
plt.legend(loc='upper right', frameon=True, framealpha=0.9)
#plt.savefig('data/rejection.svg')
```


```python

```

7c7195b48cd301a91b1100429c287ff14c3a20f1
ipynb
Jupyter Notebook
4-assets/BOOKS/Jupyter-Notebooks/Overflow/30_RejectionSampling.ipynb
impastasyndrome/Lambda-Resource-Static-Assets
7070672038620d29844991250f2476d0f1a60b0a
MIT
4-assets/BOOKS/Jupyter-Notebooks/Overflow/30_RejectionSampling.ipynb
impastasyndrome/Lambda-Resource-Static-Assets
7070672038620d29844991250f2476d0f1a60b0a
MIT
4-assets/BOOKS/Jupyter-Notebooks/Overflow/30_RejectionSampling.ipynb
impastasyndrome/Lambda-Resource-Static-Assets
7070672038620d29844991250f2476d0f1a60b0a
MIT
2021-11-05T07:48:26.000Z
2021-11-05T07:48:26.000Z
Qwen/Qwen-72B
1. YES
2. YES