Monte Carlo 在 8 个维度中未给出正确答案

Question

我一直在尝试使用蒙特卡罗方法计算 (10*6)*sin(x0+x1+x2+x3+x4+x5+x6+x7)dx0dx1dx2dx3dx4dx5dx6dx7 的积分作为我计算的一部分大学课程。经过一些尝试，我设法让代码工作，但它似乎没有给出正确的答案。我得到了 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

Monte Carlo 在 8 个维度中未给出正确答案

Monte Carlo in 8 Dimensions not giving correct answer

python

multidimensional-array

integration

montecarlo