python - 如何在 Python 中使用蒙特卡罗方法计算 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 维球体的体积吗？谢谢!

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 快速获取其长度并查看该点是否位于以原点为中心的半径为 1 的球体内。

这是一些更正的代码。我删除了一个循环——你试图建立一个正确大小的向量的那个，但我们现在用 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 的更一致的答案。 .这与 theoretical answer 非常吻合的

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

其计算结果为 2.550164039877345 .

(这是对评论的回复，用于解释返回值 np.power(2.0, dim) * (count_in_sphere / iterations 。)

此例程生成随机点，其中每个坐标在区间 [-1.0, 1.0) 中.这意味着这些点均匀分布在维度 dim 的超立方体中。 .超立方体的体积是其边的乘积。每边有长度 2.0所以我们可以用 2.0 计算超立方体的体积电源 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

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