# Draw samples from the distribution:
shape, scale = 2., 2. # mean and dispersion
s = np.random.gamma(shape, scale, 1000)
# Display the histogram of the samples, along with
# the probability density function:
import matplotlib.pyplot as plt
import scipy.special as sps
count, bins, ignored = plt.hist(s, 50, normed=True)
y = bins**(shape-1)*(np.exp(-bins/scale) /
(sps.gamma(shape)*scale**shape))
plt.plot(bins, y, linewidth=2, color='r')
plt.show()