In [2]:
import sympy
import numpy as np
import matplotlib.pyplot as plt
In [3]:
# Most distributions will be too skinny-tailed
# to use rejection sampling. A good proposal
# distribution ends up being… another Pareto
# distribution!
def pareto_inv(x_m, a):
return x_m / np.random.rand()**(1/a)
def pareto_rejection(x_m, a, n=1000):
res = []
for i in range(n):
# As long as we use a proposal with
# lower alpha, it will be fatter-tailed.
# We use xm_prop = x_m and alpha_prop = alpha – 1,
#and we can analytically derive c as
# c = alpha / (alpha – 1) # same number of searches of c
x = pareto_inv(x_m, a – 1)
# Again, we can analytically derive
# f_t(x) / f_p(x) as alpha / ((alpha – 1) x), (ratio)
# so f_t(x) / (f_p(x) c) = 1/x. (rejection probability, always remember it is less than 1)
while np.random.rand() >= 1/x:
x = pareto_inv(x_m, a-1)
res.append(x)
return np.array(res)
x_m, a = 1, 4
samples = pareto_rejection(x_m, a, 10000)
x = np.arange(1, 10, .01)
fx = a / x**(a + 1)
plt.hist(samples, density=True, bins=1000)
plt.plot(x, fx)
plt.show()

In [0]: