如何在Python中用MonteCarloMethod计算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)

这里的错误是：

^{pr2}$

知道这个算法的人能不能试着运行它，告诉我应该实现什么，得到这个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一次建立它。我还用更具描述性的名称替换了变量名，并做了一些其他更改。在

^{pr2}$

注意，仅仅1000次montecarlo迭代就给出了非常糟糕的精度。您需要更多的迭代才能得到有用的答案，所以我将重复次数改为100,000。这个例程现在比较慢，但是给出了更加一致的答案2.5。这与

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，或者{}来计算超立方体的体积。在

我们生成了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

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