如何在 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 的函数一样，它 returns 一个数组——在本例中，是一个你声明的长度的向量（在本例中 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 次蒙特卡洛迭代就给出了非常差的精度。您将需要更多的迭代才能获得有用的答案，因此我将重复次数更改为 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

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

如何在 Python 中用蒙特卡洛法计算 10 维球体的体积？

How to compute volume of 10-Dimentional sphere with Monte-Carlo-Method in Python?

python

numpy

math

montecarlo