9. Stochastic Methods: Exercise#
Name:
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
1. Integration#
When an explicit function like \(x_n = g(x_1,...,x_{n-1})\) is not available, you can use a binary function that takes 1 within and 0 outside of a region to measure its volume.
Measure the volume of an n-dimensional sphere using a binary function $\( f(x) = \left\{ \begin{array}{ccc} 1 & \mbox{if} & \sum_{i=1}^n x_i^2 \le 1 \\ 0 & \mbox{otherwise} & \end{array}\right. \)$
Measure the area or volume of a shape of your interest by sampling method.
2. Rejection Sampling#
Let us generate samples from Gamma distribution
$\( p(x; k, \theta) = \frac{1}{\Gamma(k)\theta^k}x^{k-1}e^{-\frac{x}{\theta}} \)\(
with the shape parameter \)k>0\( and the scaling parameter \)\theta\( using the gamma function \)\Gamma(k)$ available in scipy.special.
from scipy.special import gamma
from scipy.special import factorial
Here is how the gamma function looks like, together with factorial at integer values.
kmax = 5
x = np.linspace(0, kmax, 50)
plt.plot(x, gamma(x))
plt.plot(range(1,kmax+1), [factorial(k-1) for k in range(1,kmax+1)], 'o')
plt.xlabel("x"); plt.ylabel("f(x)")
Text(0, 0.5, 'f(x)')
Define the Gamma density function with arbitrary \(k\) and \(\theta\).
def p_gamma(x, k=1, theta=1):
"""density function of gamma distribution"""
return
Consider the exponential distribution \(q(x;\mu)=\frac{1}{\mu}e^{-\frac{x}{\mu}}\) as the proposal distribution.
def p_exp(x, mu=1):
"""density function of exponential distribution"""
return
def x_exp(n=1, mu=1):
"""sample from exponential distribution"""
y = np.random.random(n) # uniform in [0,1]
return
The ratio of the target and proposal distributions is $\( \frac{p(x;k,\theta)}{q(x;\mu)} = \frac{\mu}{\Gamma(k)\theta^k}x^{k-1}e^{(\frac{1}{\mu}-\frac{1}{\theta})x}. \)\( By setting \)\( \frac{d}{dx}\frac{p(x;k,\theta)}{q(x;\mu)} = 0 \)\( we have \)\( \{(k-1)x^{k-2}+(\frac{1}{\mu}-\frac{1}{\theta})x^{k-1}\}e^{(\frac{1}{\mu}-\frac{1}{\theta})x} = 0. \)\( Thus at \)\( x = \frac{(k-1)\mu\theta}{\mu-\theta} \)\( the ratio \)\frac{p(x;k,\theta)}{q(x;\mu)}\( takes the maximum \)\( \frac{\mu^k}{\Gamma(k)\theta}\left(\frac{k-1}{\mu-\theta}\right)^{k-1}e^{1-k}. \)\( 2) What is a good choice of \)\mu\( to satisfy \)p(x)\le cq(x)\( and what is the value of \)c$ for that?
By setting \(\mu=k\theta\), we have \(p(x)\le cq(x)\) with \(c=\frac{k^k}{\Gamma(k)}e^{1-k}\)
Verify that \(cq(x)\) covers \(p(x)\) by plotting them for some \(k\) and \(\theta\).
k = 2
theta = 2
c = (k**k)/gamma(k)*np.exp(1-k)
x = np.linspace(0, 10, 50)
plt.plot(x, p_gamma(x, k, theta))
plt.plot(x, c*p_exp(x, k*theta))
plt.xlabel("x");
---------------------------------------------------------------------------
ValueError Traceback (most recent call last)
Cell In[6], line 5
3 c = (k**k)/gamma(k)*np.exp(1-k)
4 x = np.linspace(0, 10, 50)
----> 5 plt.plot(x, p_gamma(x, k, theta))
6 plt.plot(x, c*p_exp(x, k*theta))
7 plt.xlabel("x")
File /opt/anaconda3/lib/python3.12/site-packages/matplotlib/pyplot.py:3590, in plot(scalex, scaley, data, *args, **kwargs)
3582 @_copy_docstring_and_deprecators(Axes.plot)
3583 def plot(
3584 *args: float | ArrayLike | str,
(...)
3588 **kwargs,
3589 ) -> list[Line2D]:
-> 3590 return gca().plot(
3591 *args,
3592 scalex=scalex,
3593 scaley=scaley,
3594 **({"data": data} if data is not None else {}),
3595 **kwargs,
3596 )
File /opt/anaconda3/lib/python3.12/site-packages/matplotlib/axes/_axes.py:1724, in Axes.plot(self, scalex, scaley, data, *args, **kwargs)
1481 """
1482 Plot y versus x as lines and/or markers.
1483
(...)
1721 (``'green'``) or hex strings (``'#008000'``).
1722 """
1723 kwargs = cbook.normalize_kwargs(kwargs, mlines.Line2D)
-> 1724 lines = [*self._get_lines(self, *args, data=data, **kwargs)]
1725 for line in lines:
1726 self.add_line(line)
File /opt/anaconda3/lib/python3.12/site-packages/matplotlib/axes/_base.py:303, in _process_plot_var_args.__call__(self, axes, data, *args, **kwargs)
301 this += args[0],
302 args = args[1:]
--> 303 yield from self._plot_args(
304 axes, this, kwargs, ambiguous_fmt_datakey=ambiguous_fmt_datakey)
File /opt/anaconda3/lib/python3.12/site-packages/matplotlib/axes/_base.py:460, in _process_plot_var_args._plot_args(self, axes, tup, kwargs, return_kwargs, ambiguous_fmt_datakey)
457 # Don't allow any None value; these would be up-converted to one
458 # element array of None which causes problems downstream.
459 if any(v is None for v in tup):
--> 460 raise ValueError("x, y, and format string must not be None")
462 kw = {}
463 for prop_name, val in zip(('linestyle', 'marker', 'color'),
464 (linestyle, marker, color)):
ValueError: x, y, and format string must not be None
Implement a function to general samples from Gamma distribution with arbitorary \(k\) and \(\theta\) using rejection sampling from exponential distribution.
def x_gamma(n=1, k=1, theta=1):
"""sample from gamma distribution by rejection sampling"""
c = (k**k)/gamma(k)*np.exp(1-k)
#print("c =", c)
xe = x_exp(n, k*theta)
paccept =
accept = np.random.random(n)<paccept
xg = xe[accept] # rejection sampling
#print("accept rate =", len(xg)/n)
return(xg)
k = 2
theta = 2
# sample histogram
xs = x_gamma(1000, k, theta)
plt.hist(xs, bins=20, density=True)
# compare with the density function
x = np.linspace(0, 10, 100)
plt.plot(x, p_gamma(x, k, theta))
plt.xlabel("")
3. Importance Sampling#
You have \(m\) sample values \(f(x_i)\) \((i=1,...m)\) at \(x_i\) following a normal distribution $\( q(x;\mu_0,\sigma_0) = \frac{1}{\sqrt{2\pi \sigma_0^2}}e^{-\frac{(x-\mu_0)^2}{2\sigma_0^2}}. \)$
Consider estimating the mean of \(f(x)\) for samples with a different normal distribution \(p(x;\mu_1, \sigma_1)\) by importance sampling $\( E_p[h(x)] = E_q\left[\frac{p(x)}{q(x)}h(x)\right] \)$
What is the importance weight \(\frac{p(x)}{q(x)}\)?
Let us consider the simplest case of \(f(x)=x\).
Generate \(m=100\) samples with \(\mu_0=100\) and \(\sigma_0=20\) and take the sample mean \(E_q[f(x)]\).
m = 100
mu0 = 100
sig0 = 20
xs =
# show histogram
plt.hist(xs, bins=20, density=True)
plt.xlabel("x")
# check the sample mean
np.mean(xs)
Estimate the mean \(E_p[f(x)]\) for \(\mu_1=120\) and \(\sigma_1=10\) by importance sampling.
mu1 = 120
sig1 = 20
importance =
mean1 = np.dot(importance, xs)/m
print(mean1)
See how the result changes with different settings of \(\mu_1\), \(\sigma_1\) and sample size \(m\).
Optional) Try with a different function \(f(x)\).
4. MCMC#
Try applying Metropolis sampling to your own (unnormalized) distribution.