import math


# gamma(n) = (n-1)!
#

print( [ math.gamma(x) for x in range(1,10) ] )

# Works for fractions too !
print( [ math.gamma(x+0.5) for x in range(1,10) ] )


# n! m!
# ----- = n choose m
# (n+m)!

def B(alpha, beta):
    """a normalizing constant so that the total probability 
       of beta_pdf is 1"""
    return math.gamma(alpha) * math.gamma(beta) / math.gamma(alpha + beta)

# https://en.wikipedia.org/wiki/Beta_distribution#Probability_density_function
#
# This is the beta distribution
#
#   f(x; a, b) = x^(a-1) (1-x)^(b-1) / B(a,b)
#
def beta_pdf(x, alpha, beta):
    if x < 0 or x > 1: # no weight outside of [0, 1]
        return 0
    return x ** (alpha - 1) * (1 - x) ** (beta - 1) / B(alpha, beta)