Решетка su(2) калибровочная теория и генерация случайных чисел в питоне

Question

В настоящее время я пишу простую программу на python для моделирования 1 + 1 мерной теории SU(2) Янга Миллса. Для случая SU(2) существует конкретный алгоритм тепловой ванны для обновления переменных Link. Однако для реализации этого алгоритма мне нужно сгенерировать случайное действительное число X, чтобы X распределялось в соответствии с P(x) = sqrt(1-X^2)*e^(k*X)где k - действительное число от отрицательной бесконечности до бесконечности.

К счастью, существует алгоритм для генерации X в соответствии с указанным распределением. Используя мои ограниченные навыки в Python, я реализовал такой алгоритм. Вот код Я использую только NumPy здесь.

def algorithm(k):
    count = 0
    while  1 != 0:
        r1,r2,r3,r4 = np.random.uniform(low=0,high=1),np.random.uniform(low=0,high=1),np.random.uniform(low=0,high=1),np.random.uniform(low=0,high=1)
        L1 = -1/(2*k)*(np.log(r1)+np.log(r3)*np.cos(2*np.pi*r2)**2)
        if r4**2 <= 1 - L1:
            X = 1 -2*L1
            break
        else:
            count = count +  1
            continue
    print(count)
    return X  

В основном, если мы берем три равномерно распределенных случайных числа в интервалах от 0 до 1, мы можем генерировать случайную величину l1, которая является функцией трех случайных чисел.

Мы принимаем это значение L1, если 1 - L1 больше или равно четвертому случайному числу в квадрате (равномерно распределенному в интервале от 0 до 1). Иначе мы возвращаемся к началу и делаем все заново. Мы делаем это, пока не примем значение L1. После того, как мы примем L1, мы вычисляем X как 1 - 2*L1. Этот алгоритм гарантирует, что X следует требуемому распределению.

В моей программе мне нужно будет сгенерировать двумерный массив X. Это довольно медленно в моей текущей реализации. Итак, вот мой вопрос; Есть ли более простой способ сделать это с помощью каких-либо предустановленных пакетов NumPy? Если такого метода не существует, есть ли способ векторизовать эту функцию для генерации двумерной решетки случайных X без простой итерации с циклом for?

Answer 1

Я не знаю, существует ли встроенная функция, которая возвращает именно то распределение, которое вы хотите, но я считаю, что векторизация вашего кода не должна быть сложной. Просто сделай r1, r2, r3 а также r4 векторы, использующие size параметр uniform функционировать как упомянуто Уорреном и выполнять эти операции. Как уже упоминалось DSM, вы также можете просто использовать кортеж как size параметр и сделать все одним вызовом.

Вы можете держать цикл и повторять операции каким-либо образом, пока не получите N значения, но я бы просто удалил цикл и оставил только те числа, которые удовлетворяют условию. Это дает меньше, чем N цифры, но это просто для кода.

Примерно так может быть то, что вы хотите:

def algorithm_2(k, N):
    r1,r2,r3,r4 = np.random.uniform(low=0,high=1, size=(4,N))
    L1 = -1/(2*k)*(np.log(r1)+np.log(r3)*np.cos(2*np.pi*r2)**2)
    reduced_L1 = L1[r4**2 <= 1 - L1]
    return 1-2*reduced_L1

Запуск это дает:

>>> algorithm_2(1, 50)
array([-0.21110547, -0.70285195,  0.0475383 , -0.20860877, -0.07776909,
       -0.21907097,  0.70566776,  0.3207524 ,  0.71130986,  0.45789795,
        0.15865965, -0.13757757,  0.04306286,  0.46003952])

Если вы хотите функцию, которая всегда возвращает точно N-векторным вектором вы можете написать оболочку, которая продолжает вызывать вышеуказанную функцию, а затем объединяет массивы. Что-то вроде этого:

def algorithm_3(k, N):
    total_size =0
    random_arrays = []
    while total_size < N:
        random_array = algorithm_2(k, N)
        total_size += len(random_array)
        random_arrays.append(random_array)
    return np.hstack(random_arrays)[:N]