from math import isnan

from scipy import stats

def beta2(x, y):
    return (gamma(x) * gamma(y)) / gamma(x + y)

vals = [beta2(a, b) for a in srange(-2, 2, 0.2) for b in srange(-2, 2, 0.2)] # test values to work out contours
filtered = filter(lambda num: (not isnan(num)) and (-1000 < num < 1000), vals)
contours = stats.mstats.mquantiles(filtered, srange(0, 1, 1/15))
contour_plot(beta2, (x, -3, 3), (y, -3, 3), cmap='rainbow', contours=contours)