In [1]:
import sympy
import numpy as np
import matplotlib.pyplot as plt
In [5]:
# We can calculate the inverse analytically:
# Integration yields F(x) = 1 – (x_m / x)**alpha.
# Solving $u= F(x)$ yields $F^{-1}(u) = x_m/(1 – u)**(1/alpha),
# which we can simplify to $F^{-1}(u) = x_m / u**(1/alpha).
#Alternatively, we can use sympy.
x = sympy.Symbol(“x”)
a = sympy.Symbol(“a”)
x_m = sympy.Symbol(“xm”)
fx = a / (x**(a+1))
Fx = sympy.integrate(fx, (x,x_m, x))
u = sympy.Symbol(‘u’)
# The given solution is odd, but parsing it we
# can see it gives the same solution.
print(sympy.solve(Fx – u, x))
# All that’s left to do is implement the algorithm!
def pareto_inv(x_m, a, n=10000):
res = []
for i in range(n):
u = np.random.rand()
res.append(x_m / u**(1/a))
return np.array(res)
x_m, a = 1, 4
samples = pareto_inv(x_m, a)
t = np.arange(x_m, 10, .01)
ft = a / t**(a + 1)
plt.hist(samples, bins=1000, density=True)
plt.plot(t, ft)
plt.show()
[Piecewise((xm*exp(u/a), Eq(a, 0)), (nan, True)), Piecewise(((-xm**a/(u*xm**a – 1))**(1/a), (a > -oo) & (a < oo) & Ne(a, 0)), (nan, True))]  In [0]: