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

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，它是三个随机数的函数。

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

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

Answer 1

我不知道是否存在一个内置函数来准确返回您想要的分布，但我相信矢量化您的代码应该不难。只需使用 Warren 提到的uniform函数的size参数制作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-ary 向量的函数，你可以编写一个包装器来不断调用上面的函数，然后连接数组。像这样：

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]