Add normal distribution in the random module
I think this issue should be splitted in two issues:
- Fast normal distrib (described below)
- Random number in any distribution (a bit more complexe... and costly, probably not what you want for the next release)
A fast method to get a random number that follows a normal distribution is to cumulate several random numbers that follows a uniform distribution.
def random_normal(mu, sigma):
N = 10 # Number of uniform r.v. to cumulate
r = 0
for _ in range(N):
r += np.random.random()*2-1 # Cumulate uniform r.v. in [-1,1]
return mu + sigma * r / np.sqrt(N/3)
With N> 10 or above, the obtained distributions are very satisfying:
To generate a random number in any distribution, I suggest to use a MCMC method such as the well known Metropolis algorithm.
The idea consist in placing a walker randomly in the desired area in wich we want to generate a random number. We then generate randomly a new position for the walker in this area. If the distribution is lower at the new position than at the old one, we accept the move only if a random number between 0 and 1 is lower than the ratio of the distribution at these two positions.
def metropolis(dist, a, b):
x = np.random.random() * (b - a) + a
N = 10
for _ in range(N):
x2 = np.random.random() * (b - a) + a
if np.random.random() < dist(x2) / dist(x):
x = x2
return x