Sample transformation power.py

From Werner KRAUTH

import random, math
import matplotlib.pyplot as plt

N_trials = 1000000
data = []
gamma = -0.7
for iter in range(N_trials):
    Upsilon = random.uniform(0.0, 1.0)
#
#   This is the sample transformation SMAC eqs (1.28), (1.29)
#
    x = Upsilon ** (1.0 / (gamma + 1))
    data.append(x)

plt.title('power-law distribution (sample transformation)  $\gamma = $ '+ str(gamma))
plt.xlabel('$x$')
plt.ylabel('$\pi(x)$')
plt.hist(data, bins=100, density=True,label='data')
XValues = []
YValues = []
for i in range(5, 1000):
    x = i / 1000.0
    XValues.append(x)
    YValues.append((gamma + 1.0) * x ** gamma)
plt.plot(XValues, YValues, label='theory')
plt.legend(loc='upper right')
plt.show()