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