from scipy.stats import hypergeom

n = 1000
k = 20
x = 1
lhx = 0.0
be = 0
for b in range(x, n - k + x):
    lh = hypergeom.pmf(x, n, b, k)
    if lh > lhx:
        be = b
        lhx = lh
        
print("be =",be)