0

I'm currently writing a simple program in python to simulate 1 + 1 dimensional SU(2) yang mills theory. For the case of SU(2) there exists a particular heatbath algorithm for updating the Link variables. In order to implement this algorithm however I need to generate a random real number X such that X is distributed according to P(x) = sqrt(1-X^2)*e^(k*X), where k is a real number from negative infinity to infinity.

Fortunately there exists an algorithm to generate X's according to said distribution. Using my limited skills in python I implemented such an algorithm. Here is the code. I'm only using numpy here. 

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

Basically, if we take three uniformly distributed random numbers in the intervals 0 to 1 we can generate a random variable l1 which is a function of the three random numbers.

We accept this value L1 if 1 - L1 is greater than or equal to a fourth random number squared (uniformly distributed in the interval 0 to 1). Else we loop back to the beginning and do it all over again. We do this until we accept a value of L1. After we accept L1 we compute X as being 1 - 2*L1. This algorithm ensures that X follows the required distribution.

In my program I'm going to have to generate an two dimensional array of X's. This is quite slow in my current implementation. So here's my question; is there a simpler way to do this using any preset numpy packages? If such a method doesn't exist, is there a way to vectorize this function to generate a two dimensional lattice of random X's without simply iterating it with a for loop?

4

1 回答 1

1

我不知道是否存在可以准确返回您想要的分布的内置函数，但我相信矢量化您的代码应该不难。只需使用Warren 提到的函数参数r1制作、r2和向量r3并执行这些操作。正如 DSM 所提到的，您也可以只使用元组作为参数，并通过一次调用完成所有操作。r4sizeuniformsize

您可以保留循环并以某种方式重复操作，直到获得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]
于 2018-08-06T20:55:35.473 回答