Ctrl+K

Probability distributions Exercise 1

Probability distributions Exercise 1#

import numpy as np
import matplotlib.pyplot as plt

seed = 5885
rng = np.random.default_rng(seed)

# We'll need this a few times, so write a function to take the x values and 
# plot a histogram comparing to the analytic function func
def plot_distribution(x, func):
    plt.clf()
    plt.hist(x, density=True, bins=100, histtype = 'step')
    xx = np.linspace(0.0,max(x),100)
    plt.plot(xx, func(xx),':')
    plt.yscale('log')
    plt.xlabel('x')
    plt.show()

Rejection method#

def f(x):
    return np.exp(-x)

xmax = 10
x = xmax * rng.uniform(size = 10**5)    # Generate x values up to xmax
y = rng.uniform(size = 10**5)

# Reject
ind = np.where(y <= f(x))
ind2 = np.where(y > f(x))

# Plot the points to show which ones are accepted and which rejected
plt.plot(x[ind2],y[ind2], 'bo', ms=0.1)
plt.plot(x[ind],y[ind], 'ro', ms=0.1)
# Plot the analytic solution for the boundary (see below)
xx = np.linspace(0.0,xmax,100)
plt.plot(xx, f(xx), 'k')
plt.xlabel('x')
plt.show()

plot_distribution(x[ind], f)
_images/30e55a2cfa910cdf4569ffdcc7f32b23af44de4f3ff4056c9710ef3853707538.png _images/94aae17ce5b5f6bd7fab747830e225b1d86c244502747538f73e5236da8815f7.png

Transformation method#

y = rng.uniform(size = 10**5)

x = - np.log(y)

# Plot the distribution of x and compare with the exponential distribution
plot_distribution(x, f)

# Plot 1-y(x) to illustrate that y is related to the cumulative distribution function
plt.clf()
plt.plot(x, 1-y, '.', ms=0.1)
xx = np.linspace(0.0,xmax,100)
plt.plot(xx, 1-f(xx), 'k:')
plt.xlabel('x')
plt.ylabel('1-y')
plt.show()
_images/87b6e2b364d84e4e49d217a4c3392cc5c503e5d227253d6d5dd064be35ed33e2.png _images/cd527dc37a9e371a17d5859283b6ee40cbea939f0290d42605296bae8bbec467.png

Ratio of uniforms#

def f(x):
    return np.exp(-x)

u = rng.uniform(size = 10**5)
v = rng.uniform(size = 10**5)
# choose which points to keep -- we renormalize f here to drop the factor of 2
ind = np.where(u <= np.sqrt(f(v/u)))
ind2 = np.where(u > np.sqrt(f(v/u)))
x = v[ind]/u[ind]

# Plot the points to show which ones are accepted and which rejected
plt.plot(u[ind2],v[ind2], 'bo', ms=0.1)
plt.plot(u[ind],v[ind], 'ro', ms=0.1)
# Plot the analytic solution for the boundary (see below)
uu = np.linspace(0.001,1.0,100)
plt.plot(uu, -2*uu*np.log(uu), 'k')
plt.xlabel('u')
plt.ylabel('v')
plt.show()

# Plot the distribution f(x)
plot_distribution(x, f)
_images/7a3457ed056b9a89496a570df26ce79fcc9992808ef691ab5870726e4136bb3c.png _images/70cf92cc2a10afbe54432133eabae6a02a41b1b2823592ff48cc756999275efb.png

The black curve above separates the accepted and rejected points. To derive this, start with

\[u \leq \sqrt{\exp\left(-{v\over u}\right)}\]

and invert this to get \(v\):

\[u^2 \leq \exp\left(-{v\over u}\right)\]
\[2\ln u \leq -{v\over u}\]
\[2 u \ln u \leq -v\]
\[v \leq -2 u \ln u\]

which is plotted as the black curve in the figure.

Contents