如何加速python中的随机数组生成

Question

我想使用Monte-Carlo积分方法找到n球的体积并将其与分析结果进行比较.我想这样做10 ^ 6点和10 ^ 9点,虽然它有点10 ^ 6点(n = 2(圆圈),n = 3(球体)和n = 12需要大约一分钟), 10 ^ 9分非常慢.

MC方法的简短说明:为了找到半径r = 1的n球的体积,我想象一个简单的已知体积(比如一个边长为2*r的n立方体),它完全包含n球.然后我从n立方体中的均匀分布点进行采样,并检查该点是否位于球体中.我统计了n球内所有生成的点.V_sphere/V_cube的比率可以近似为N_inside/N_total,因此,V_sphere = V_cube*N_inside/N_total

这是功能:

def hyp_sphere_mc(d,samples):    

    inside = 0                     #number of points inside sphere                                             
    sum = 0                        #sum of squared components

    for j in range(0,samples):        

        x2 = np.random.uniform(0,1,d)**2     #generate random point in d-dimensions                           
        sum = np.sum(x2)                     #sum its components

        if(sum < 1.0):                                              
            inside += 1                      #count points inside sphere

    V = ((2)**d)*(float(inside)/samples)     #V = V_cube * N_inside/N_total           

    V_true = float(math.pi**(float(d)/2))/math.gamma(float(d)/2 + 1) #analytical result 

    ERR = (float(abs(V_true-V))/V_true)*100        #relative Error

    print "Numerical:", V, "\t" , "Exact: ", V_true, "\t", "Error: ", ERR

我想问题是,对于每次迭代,我生成一个新的随机数组,这需要花费很多时间,特别是如果我有10 ^ 9次迭代.有什么方法可以加快速度吗？

Answer 1

您可以使用以下内容替换循环:

inside = np.sum(np.sum(np.random.rand(samples,d)**2,1)<1)

当使用numpy时,你应该尽量避免循环.我们的想法是,您可以在矩阵中一次生成所有样本,然后对所有后续操作进行矢量化.