import matplotlib.pyplot as plt
import numpy as np
from scipy.special import betaln, gammaln
def factln(x):
    return gammaln(x + 1)
def logp(x, alpha, beta, r):
    return (
        betaln(r + x, alpha + beta)
        - betaln(r, alpha)
        + gammaln(x + beta)
        - factln(x)
        - gammaln(beta)
    )
plt.style.use('arviz-darkgrid')
x = np.arange(0, 25)
params = [
    (1, 1, 1),
    (1, 1, 10),
    (1, 10, 1),
    (1, 10, 10),
    (10, 10, 10),
]
for alpha, beta, r in params:
    pmf = np.exp(logp(x, alpha, beta, r))
    plt.plot(x, pmf, "-o", label=r'$alpha$ = {}, $beta$ = {}, $r$ = {}'.format(alpha, beta, r))
plt.xlabel('x', fontsize=12)
plt.ylabel('f(x)', fontsize=12)
plt.legend(loc=1)
plt.show()