Example 2.13#
As stated in notes, we aim at comparing two procedures. First,
Sample
,Sample
and comparing resulting
.
By analytical results, they must match.
import numpy as np
import matplotlib.pyplot as plt
rng = np.random.default_rng(12)
mu = 2
sigma = 0.5
N = 100000
x = np.zeros(N)
y = np.zeros(N)
# We use rng.normal for this
# Note that rng.normal takes the standard deviation as the second argument
# You can use Box-M\"uller transform to sample Gaussians too
for i in range(N):
x[i] = rng.normal(mu, sigma, 1)
y[i] = rng.normal(x[i], 1, 1)
def marginal_y(z, mu, sigma):
return np.exp(-0.5 * (z-mu)**2/(sigma**2 + 1)) / np.sqrt(2 * np.pi * (sigma**2 + 1))
yy = np.linspace(-2, 7, 1000)
plt.hist(y, bins=100, density=True, label='y', color='black')
plt.plot(yy, marginal_y(yy, mu, sigma), label='p(y)', color=[0.8, 0, 0])
plt.legend()
plt.show()
<ipython-input-1-e73d65d55468>:18: DeprecationWarning: Conversion of an array with ndim > 0 to a scalar is deprecated, and will error in future. Ensure you extract a single element from your array before performing this operation. (Deprecated NumPy 1.25.)
x[i] = rng.normal(mu, sigma, 1)
<ipython-input-1-e73d65d55468>:19: DeprecationWarning: Conversion of an array with ndim > 0 to a scalar is deprecated, and will error in future. Ensure you extract a single element from your array before performing this operation. (Deprecated NumPy 1.25.)
y[i] = rng.normal(x[i], 1, 1)