Boilerplate code for working with a Gamma Distribution using scipy that matches the functional form defined on Wikipedia: https://en.wikipedia.org/wiki/Gamma_distribution

gamma_example.py

import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns  # optional, just for aesthetics
from scipy.stats import gamma

# recreate the example shown in Wikipedia for reference
# param = (alpha, beta) == (k, 1/theta)
params = [
    (1, 1/2),
    (2, 1/2),
    (3, 1/2),
    (5, 1),
    (9, 2),
    (7.5 ,1),
    (0.5, 1)
]

fig, ax = plt.subplots(figsize=(8, 5))

for alpha, beta in params:
    distribution = gamma(a=alpha, scale=1/beta)
    
    x = np.linspace(0, 20, num=1000)
    pdf = distribution.pdf(x)
    
    label = f'$\\alpha={alpha}, \\beta={beta}$'
    ax.plot(x, pdf, linewidth=2, label=label)
    
    
ax.set_ylim(0, 0.5)
ax.set_xlabel('$x$')
ax.set_ylabel('$pdf(x)$')
ax.legend(frameon=False)

# optional for aesthetics
sns.despine(trim=True, offset=10)

plt.show()