How to generate random numbers from a gamma distribution with python ?

Generating random numbers from a gamma distribution with python can be done through several methods:

Using numpy.random.gamma()

The most popular way is to use the NumPy library and its built-in function, numpy.random.gamma(). This function takes in two parameters - shape & scale - which are used to define the shape of the Gamma distribution that will be generated. Additionally, you can also set a size parameter to generate multiple random numbers at once.

import numpy as np

shape = 3.0
scale = 1.0

nb_random_points = 10000

data = np.random.gamma(shape, scale, nb_random_points)

Output

array([3.14118269, 5.10012155, 1.90110229, ..., 1.66629833, 0.51727531,
   0.37046791])



import matplotlib.pyplot as plt
import numpy as np

plt.hist(data, 100, density=True)

plt.title('How to generate random numbers \n from a gamma distribution with python ?')

plt.savefig("gamma_distribution_random_numbers_01.png", bbox_inches='tight',dpi=200, facecolor='white')

plt.show()

Simulate a Gamma distribution using scipy’s stats module

Alternatively, you can simulate a Gamma distribution with Python using scipy’s stats module and its method called gamma(). This method accepts three parameters - shape, scale, and size - and it works similarly as the numpy method mentioned earlier.

import matplotlib.pyplot as plt
import numpy as np
import scipy.special as sps

plt.hist(data, 100, density=True)

X = np.linspace(0,16,100)

PDF = ( 1.0 / ( sps.gamma(shape) * scale**shape) ) * X**(shape-1) * np.exp(-X/scale)

plt.plot(X,PDF)

plt.title('How to generate random numbers \n from a gamma distribution with python ?')

plt.savefig("gamma_distribution_random_numbers_02.png", bbox_inches='tight',dpi=200, facecolor='white')

plt.show()

Generate random numbers using the cumulative distribution function

To calculate the cumulative distribution function (CDF) of a gamma function a solution is to use scipy stats:

from scipy.stats import gamma

plt.hist(data, 100, density=True)

X = np.linspace(0,16,1000)

CDF = gamma.cdf(X, shape, 0, scale)

plt.plot(X,CDF)

plt.title('How to generate random numbers \n from a gamma distribution with python ?')

plt.xlim(0,16)

plt.savefig("gamma_distribution_random_numbers_03.png", bbox_inches='tight',dpi=200, facecolor='white')

plt.show()

Generating random numbers using the cumulative distribution function (CDF) involves assigning a probability to each outcome and then using the CDF to determine which value is chosen. To generate such a number, a user must first define the probability of each outcome, typically as a range from 0 to 1:

from scipy.optimize import minimize

x = np.random.random_sample((10000,))

y = np.zeros(x.shape)

def diff(x,a):
    yt = gamma.cdf(x, shape, 0, scale)
    return (yt - a )**2

for idx,x_value in enumerate(x):
    res = minimize(diff, 2, args=(x_value), method='Nelder-Mead', tol=1e-6)
    y[idx] = res.x[0]

Now let's create an histogram of random numbers using matplotlib:

plt.hist(y, 100, density=True)

X = np.linspace(0,16,1000)

X = np.linspace(0,16,100)

PDF = ( 1.0 / ( sps.gamma(shape) * scale**shape) ) * X**(shape-1) * np.exp(-X/scale)

plt.plot(X,PDF)

plt.savefig("gamma_distribution_random_numbers_04.png", bbox_inches='tight',dpi=200, facecolor='white')

plt.show()

References