scipy.special.beta

scipy.special.beta(x1, x2[, out])

y=beta(a,b) returns gamma(a) * gamma(b) / gamma(a+b)