8维蒙特卡罗没有给出正确答案

Question

我一直在用蒙特卡罗方法计算（10*6）*sin（x0+x1+x2+x3+x4+x5+x6+x7）dx0dx1dx2dx3dx4dx5dx6dx7的积分，这是我在uni的计算课程的一部分。经过一番尝试，我成功地让代码运行起来，但它似乎没有给出正确的答案。我得到了一个参考分析结果537.187。。。。但我的代码收敛到大约360

我不知道还有什么办法来解决这个问题。这与我的随机数有关吗？它们可以通过单独定义来关联吗？我应该使用单个8维生成器吗？还是我需要使用更分层的取样器

抱歉，如果是小问题，我对Python还不太熟悉

更新：很抱歉，我忘了说我使用了什么限制。我在每个维度上从0积分到pi/8

a=0          #Defining Limits of Integration
b=(np.pi)/8
N=20000    #Number of random numbers generated in each dimension

x0rand, x1rand, x2rand, x3rand, x4rand, x5rand, x6rand, x7rand = [[0]*N]*8
for i in range(N):
   x0rand[i]= np.random.uniform(a,b) #Generates N random numbers 8D between 0 and pi/8
   x1rand[i]= np.random.uniform(a,b)
   x2rand[i]= np.random.uniform(a,b)
   x3rand[i]= np.random.uniform(a,b)
   x4rand[i]= np.random.uniform(a,b)
   x5rand[i]= np.random.uniform(a,b)
   x6rand[i]= np.random.uniform(a,b)
   x7rand[i]= np.random.uniform(a,b)

def func(x0,x1,x2,x3,x4,x5,x6,x7):  #Defining the integrand
return (10**6)*np.sin(x0+x1+x2+x3+x4+x5+x6+x7)

integral = 0.0
for i in range(N):
   integral += func(x0rand[i],x1rand[i],x2rand[i],x3rand[i],x4rand[i],x5rand[i],x6rand[i],x7rand[i])

answer=(((b-a)**8)/N)*integral
print ('The Integral is:', answer)

Answer 1

这些代码收敛于您的分析答案。我通过一次生成所有样本使其更加紧凑。基本上，您生成一个随机数矩阵，每行是一个样本

a=0          #Defining Limits of Integration
b=(np.pi)/8
N=200000    #Number of random numbers generated in each dimension

samples = np.random.uniform(a,b,(N,8))

def func(x):  #Defining the integrand
    return (10**6)*np.sin(np.sum(x))

integral = 0.0
for i in range(N):
    integral += func(samples[i,:])

answer=(((b-a)**8)/N)*integral
print ('The Integral is:', answer)

编辑： 更快更高效的版本

a=0          #Defining Limits of Integration
b=(np.pi)/8
N=2000000    #Number of random numbers generated in each dimension

samples = np.random.uniform(a,b,(N,8))

integrand_all_samples = (10**6)*np.sin(np.sum(samples, axis=1))
sum_all_samples = np.sum(integrand_all_samples)


answer=(((b-a)**8)/N)*sum_all_samples
print ('The Integral is:', answer)

编辑2：出了什么问题

有问题的代码可以减少为：

N=4

a,b = [[0]*N]*2
# [[0]*N]*2: [[0,0,0,0],[0,0,0,0]]
# a: [0,0,0,0]
# b: [0,0,0,0]
for i in range(N):
    a[i]= 1
    b[i] = 2
# a: [2,2,2,2]
# b: [2,2,2,2]

a和b指向内存中的同一列表，尽管名称不同。您可以通过查看id（a）和id（b）的输出来检查这一点。你真的只有一张单子。您可以在此处找到原因的详细信息：List of lists changes reflected across sublists unexpectedly