如何用Python中的Monte-Carlo方法计算10维球体积？

Question

我试图用python计算10维球体的体积,但我的计算不起作用.

这是我的代码:

def nSphereVolume(dim,iter):
    i = 0
    j = 0
    r = 0
    total0 = 0
    total1 = 0

    while (i < iter):
        total0 = 0;
        while (j < dim):
            r = 2.0*np.random.uniform(0,1,iter) - 1.0

            total0 += r*r
            j += 1

        if (total0 < 1):
            total1 += 1
        i += 1

    return np.pow(2.0,dim)*(total1/iter)

nSphereVolume(10,1000)

这里的错误:

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-253-44104dc9a692> in <module>()
     20     return np.pow(2.0,dim)*(total1/iter)
     21 
---> 22 nSphereVolume(10,1000)

<ipython-input-253-44104dc9a692> in nSphereVolume(dim, iter)
     14             j += 1
     15 
---> 16         if (total0 < 1):
     17             total1 += 1
     18         i += 1

ValueError: The truth value of an array with more than one element is ambiguous. Use a.any() or a.all()

知道这个算法的人是否可以尝试运行它并告诉我实现什么是正确的,以获得这个10-Dim球体的体积？谢谢!

Answer 1

您的日常工作中有多个问题。

您收到的错误消息来自您的行

r = 2.0*np.random.uniform(0,1,iter) - 1.0

函数调用np.random.uniform(0,1,iter)不会创建单个随机数。而是像大多数numpy函数一样，它返回一个数组-在这种情况下，是您声明的长度的向量（在这种情况下iter）。所以r是一个数组为好，因为它使用这个数组，total0是出于同样的原因数组。

然后，稍后您尝试评估total0 < 1。但是左边是一个数组，右边是一个标量，所以numpy不喜欢比较。我可以更详细地说明错误消息的含义，但这是基本思想。

解决方案是将其r视为向量，实际上是将其作为所需边长的球体中的随机点2。您输入了错误并使用了错别字，iter而不是dim将其用作随机向量的大小，但是如果进行更改，您将得到想要的点。您可以使用numpy快速获取其长度，并查看该点是否位于以原点为中心的半径球体内。

这是一些更正的代码。我删除了一个循环-在该循环中，您尝试构建正确大小的向量，但是现在numpy一次构建了所有循环。我还用更具描述性的名称替换了您的变量名称，并进行了其他一些更改。

import numpy as np

def nSphereVolume(dim, iterations):
    count_in_sphere = 0

    for count_loops in range(iterations):
        point = np.random.uniform(-1.0, 1.0, dim)
        distance = np.linalg.norm(point)
        if distance < 1.0:
            count_in_sphere += 1

    return np.power(2.0, dim) * (count_in_sphere / iterations)

print(nSphereVolume(10, 100000))

请注意，仅进行1000次Monte-Carlo迭代就会产生非常差的精度。您将需要进行更多次迭代才能获得有用的答案，因此我将重复次数更改为100,000。现在，该例程速度较慢，但​​给出了大约一致的答案2.5。这与一致以及理论答案的

np.pi**(10 // 2) / math.factorial(10 // 2)

计算结果为2.550164039877345。

（这是对评论的回应，以解释返回的值np.power(2.0, dim) * (count_in_sphere / iterations。）

此例程生成随机点，其中每个坐标都在interval中[-1.0, 1.0)。这意味着这些点在尺寸的超立方体中均匀分布dim。超立方体的体积是其侧面的乘积。每一边都有长度，2.0因此我们可以使用2.0power dim或计算超立方体的体积np.power(2.0, dim)。

我们生成了iterations点，其中count_in_sphere的点1位于以原点为中心的半径超球面中。因此，超立方体中也属于超球体的点的分数或比例为count_in_sphere / iterations。这些点在超立方体中均匀分布，因此我们可以估计，超球体体积与超立方体体积的分数与这些集合中随机点的分数相同。因此，按高中比例，我们得到

volume_of_hypersphere / volume_of_hypercube = points_in_hypersphere / points_in_hypercube

意识到方程只是近似的。将该等式的两边乘以volume_of_hypercube，我们得到

volume_of_hypersphere = volume_of_hypercube * points_in_hypersphere / points_in_hypercube

代替，我们得到

volume_of_hypersphere = np.power(2.0, dim) * count_in_sphere / iterations

这是函数返回的值。同样，它只是近似值，但是概率论告诉我们，我们生成的随机点越多，近似值越好。