Использование векторизации для генерации правильно расположенных случайных массивов

Question

Я работаю над классом статистической механики в python через Coursera. Я пытался решить проблему, когда мы видим, сколько конкретных допустимых состояний занято. Я чувствую, что мой метод должен быть математически эквивалентен приведенному, но он не дает мне одинаковых ответов (статистически, число занятых состояний слишком мало в моем методе), и он идет медленнее. Может кто-нибудь помочь мне с тем, что мне не хватает?

Две функции:

def disk_sampling_theirs(N, sigma):
    condition = False
    while condition == False:
        L = [(random.uniform(sigma, 1.0 - sigma), random.uniform(sigma, 1.0 - sigma))]
        # generate a random pair of points
        for k in range(1, N):
            #generate a new second pair of points
            a = (random.uniform(sigma, 1.0 - sigma), random.uniform(sigma, 1.0 - sigma))

            #check to see that none of the points are within 2*sigma of any of the pairs in L
            min_dist = min(math.sqrt((a[0] - b[0]) ** 2 + (a[1] - b[1]) ** 2) for b in L)

            if min_dist < 2.0 * sigma:
                #if the points are too close, throw everything out and start over with new L
                condition = False
                break
            else:
                # add a to L, and go to next point.  If condition = True (stays true) for the whole set of points, then you can end the while loop.
                L.append(a)
                condition = True
    return L

Моя векторизованная версия с использованием linalg.norm:

def disk_sampling_mine(N, sigma):
    condition = False
    while condition == False:
        L = np.random.uniform(low=sigma, high=1-sigma, size=(1,2))
        for _ in range(1,N):
            a = np.random.uniform(low=sigma, high=1-sigma, size=(1,2))
            L = np.append(L,a, axis=0)
            min_dist = min((np.linalg.norm(p[0] - p[1]) 
                           for p in itertools.combinations(L, r=2)))
            if min_dist < 4 * sigma**2:
                condition = False
                break
            else:
                condition = True
    return list(L)

Затем я проверяю код, используя скрипт, который они дали

sigma = 0.15
del_xy = 0.05
n_runs = 100000
conf_a = ((0.30, 0.30), (0.30, 0.70), (0.70, 0.30), (0.70,0.70))
conf_b = ((0.20, 0.20), (0.20, 0.80), (0.75, 0.25), (0.75,0.75))
conf_c = ((0.30, 0.20), (0.30, 0.80), (0.70, 0.20), (0.70,0.70))
configurations = [conf_a, conf_b, conf_c]
hits = {conf_a: 0, conf_b: 0, conf_c: 0}
for run in range(n_runs):
    x_vec = disk_sampling_theirs(4, sigma)
    for conf in configurations:
        condition_hit = True
        for b in conf:
            condition_b = min(max(abs(a[0] - b[0]), abs(a[1] - b[1])) for a in x_vec) < del_xy
            condition_hit *= condition_b
        if condition_hit:
            hits[conf] += 1

for conf in configurations:
    print (conf, hits[conf])

И для моего трио состояний я получаю: (12, 20, 17), их метод (примерно, они немного меняются) (1 , 0, 0) мой метод (я редко получаю более 1 или 2 хитов)

Чего мне не хватает?