import matplotlib.pyplot as plt
import numpy as np
import scipy.stats as st
import arviz as az
plt.style.use('arviz-darkgrid')
x = np.arange(-15, 15)
params = [
    (1, 1),
    (5, 5),
    (5, 1),
]
for mu1, mu2 in params:
    pmf = st.skellam.pmf(x, mu1, mu2)
    plt.plot(x, pmf, "-o", label=r'$\mu_1$ = {}, $\mu_2$ = {}'.format(mu1, mu2))
plt.xlabel('x', fontsize=12)
plt.ylabel('f(x)', fontsize=12)
plt.legend(loc=1)
plt.show()