我正在尝试使用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维球体的体积吗?谢谢!

最佳答案

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

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

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 power dimnp.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 - 如何在Python中使用蒙特卡罗方法计算10维球体的体积?,我们在Stack Overflow上找到一个类似的问题:https://stackoverflow.com/questions/51091377/

10-12 16:34