蒙特卡洛在8个维度中没有给出正确的答案 [英] Monte Carlo in 8 Dimensions not giving correct answer

问题描述

我一直在尝试使用Monte-Carlo方法计算(10 * 6)* sin(x0 + x1 + x2 + x3 + x4 + x5 + x6 + x7)dx0dx1dx2dx3dx4dx5dx6dx7的积分，这是我的计算课程的一部分单经过一些尝试，我设法使代码正常工作，但是它似乎没有给出正确的答案.我得到的参考分析结果为537.187 ....，但是我的代码收敛到了大约360

I've been trying to compute the integral of (10*6)*sin(x0+x1+x2+x3+x4+x5+x6+x7)dx0dx1dx2dx3dx4dx5dx6dx7 using Monte-Carlo methods as part of my computing coursework at uni. I've managed to get the code to work after some trying, but it doesn't seem to give the correct answer. I've been given a reference analytical result of 537.187.... but my code converges to about 360

我不确定还可以尝试解决什么问题.这与我的随机数有关吗，是否可以通过分别定义它们进行关联?我应该使用一个8维生成器吗?还是我需要使用更分层的采样器?

I'm not sure what else to try to resolve it. Is it something to do with my random numbers, could they be correlated by defining them seperately? Should I use a single 8 dimensional generator? Or do I need to use a more stratified sampler?

很抱歉，如果它很小，我对Python还是很陌生！

Apologies if it's something small, I'm quite new to Python!

更新:对不起，我忘了说我在使用什么限制.我在每个维度上都从0集成到pi/8

UPDATE: I'm sorry I forgot to say what limits I was using. I'm integrating from 0 to pi/8 in each dimension

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)

推荐答案

此代码收敛到您的分析性答案.通过一次性生成所有样本，我使其变得更加紧凑.基本上，您会生成一个随机数矩阵，并且每一行都是一个样本.

This codes converges to your analytical answer. I've made it more compact by generating all the samples in one go. Basically, you generate a matrix of random numbers and each row is one sample.

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)

出了什么问题

有问题的代码可以简化为:

The problematic code can be reduced to:

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)的输出进行检查.您确实只有一个清单.您可以在此处找到原因的详细信息:在子列表中反映的列表更改列表意外地

a and b are pointing to the same list in memory despite having the different names. You can check this by looking at the output of id(a) and id(b). You truly only have one list. You can find the details of why here: List of lists changes reflected across sublists unexpectedly

这篇关于蒙特卡洛在8个维度中没有给出正确的答案的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！