python - 格子 su(2) 规范理论和 python 中的随机数生成

Question

我目前正在用 python 编写一个简单的程序来模拟 1 + 1 维 SU(2) yang mills 理论。对于 SU(2) 的情况，存在用于更新链接变量的特定热浴算法。但是，为了实现此算法，我需要生成一个随机实数 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，它是这三个随机数的函数。

如果 1 - L1 大于或等于第四个随机数的平方(均匀分布在 0 到 1 区间)，我们接受这个值 L1。否则我们回到开头并重新做一遍。我们这样做，直到我们接受 L1 的值。在我们接受 L1 之后，我们将 X 计算为 1 - 2*L1。该算法确保 X 服从所需的分布。

在我的程序中，我将不得不生成一个 X 的二维数组。这在我当前的实现中非常慢。所以这是我的问题；有没有更简单的方法使用任何预设的 numpy 包来做到这一点？如果不存在这样的方法，是否有一种方法可以将该函数矢量化以生成随机 X 的二维点阵，而无需简单地使用 for 循环对其进行迭代？

Answer 1

我不知道是否存在可以准确返回您想要的分布的内置函数，但我相信对您的代码进行矢量化应该不难。只需使用 size parameter of the uniform function 生成 r1、r2、r3 和 r4 向量正如沃伦所提到的那样并执行这些操作。正如 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]