如何用Python中的Monte-Carlo方法计算10维球体积？

我试图用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-Dim球体的体积？谢谢！

解决方法:

你的日常生活中有很多问题.

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

r = 2.0*np.random.uniform(0,1,iter) - 1.0

函数调用np.random.uniform(0,1,iter)不会创建单个随机数.相反,像大多数numpy函数一样,它返回一个数组 – 在这种情况下,是一个你声明的长度的向量(在本例中是iter).所以r也是一个数组,因为它使用这个数组,而total0是一个数组,原因相同.

然后你试着评估total0<但是左边是一个数组,右边是一个标量,所以numpy不喜欢比较.我可以详细了解错误消息的含义,但这是基本的想法. 解决方案是将r视为向量 – 实际上,作为球体中的随机点,您需要边长2.你做了一个拼写错误并使用iter而不是随机向量的大小变暗,但是如果你进行了改变,你将得到你想要的点.您可以使用numpy快速获取其长度,并查看该点是否位于以原点为中心的半径范围内. 这是一些纠正的代码.我删除了一个循环 – 你试图建立一个正确大小的向量的循环,但我们现在有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功率调暗或者np.power(2.0,暗淡)来计算超立方体的体积.

我们生成了迭代点,它们的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,math,numpy,montecarlo
来源： https://codeday.me/bug/20190522/1152827.html