1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
|
import numpy as np
from scipy.stats.sampling import NumericalInverseHermite
from scipy.stats import norm
from scipy.special import ndtr
import matplotlib.pyplot as plt
class StandardNormal:
def pdf(self, x):
return 1/np.sqrt(2*np.pi) * np.exp(-x**2 / 2)
def cdf(self, x):
return ndtr(x)
dist = StandardNormal()
urng = np.random.default_rng()
rng = NumericalInverseHermite(dist, random_state=urng)
rvs = rng.rvs(10000)
x = np.linspace(rvs.min()-0.1, rvs.max()+0.1, 1000)
fx = norm.pdf(x)
plt.plot(x, fx, 'r-', lw=2, label='true distribution')
plt.hist(rvs, bins=20, density=True, alpha=0.8, label='random variates')
plt.xlabel('x')
plt.ylabel('PDF(x)')
plt.title('Numerical Inverse Hermite Samples')
plt.legend()
plt.show()
|
