1

This is my first time here. Apologies, I'm quite new to python. I'm trying to use a monte carlo method to compute an integral in 8D as part of my uni coursework. However, when I run the program, I get the error described in the title on line when i define x0[i]. It seems to only allow i to run up to 2, but I want to investigate the error as N is varied so I need to be able to use N>2. Any help would be appreciated

Jake

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

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

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)

IndexError                                Traceback (most recent call last)
<ipython-input-141-1a3483405a4d> in <module>()
  4 
  5 for i in range(N):
  ----> 6   x0rand[i]=((np.pi)/8)*np.random.uniform(0,1) #Generates N random numbers 8D between 0 and pi/8
  7   x1rand[i]=((np.pi/8))*np.random.uniform(0,1)
  8   x2rand[i]=((np.pi/8))*np.random.uniform(0,1)

 IndexError: index 2 is out of bounds for axis 0 with size 2
Jake A
  • 25
  • 4

2 Answers2

0

The posted program isn't the complete code which has run, otherwise you'd have got the message NameError: name 'x0rand' is not defined. You must have defined and dimensioned x0rand before, and that with size 2 rather than N.

To (re-)define x0rand etc. appropriately, you can put

x0rand, x1rand, x2rand, x3rand, x4rand, x5rand, x6rand, x7rand = np.zeros((8, N))

before the loop.

Armali
  • 18,255
  • 14
  • 57
  • 171
  • 1
    Don't want to downvote but that answer generates quite an unexpected behaviour for a new python programmer. Especially, considering that afterwards comes a loop that will mutate its elements. More specifically, x0rand, x1rand... are all the same list! So when you mutate one you mutate all of them. This "unexpected" behaviour lead down the road to the following question from the OP (https://stackoverflow.com/questions/66704954/monte-carlo-in-8-dimensions-not-giving-correct-answer/66708811#66708811) Please consider adding some caveats to your answer. – GuillemB Mar 19 '21 at 14:14
  • Thank you for your comment; you are of course absolutely right. But hey, because of my error the querist learned something new, and I was reminded of something I should have thought of from the beginning. ;-) Rather than _adding some caveats_ to my erroneous answer, I corrected it. Thanks again! – Armali Mar 19 '21 at 19:58
  • 1
    No bother at all. Correcting the answer is clearly the better way to go. Cheers! – GuillemB Mar 20 '21 at 08:48
0

Unfortunately this code doesn't work, anyway i suggest to use the power of numpy:

Some HIT

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

can be initialized in vectorial mode, you can see size parameter.

xrand = np.random.uniform(0,1,size=(N,))

The function, itself can become more efficient. Calling it will be more readable.

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

func(xrand)
Glauco
  • 1,385
  • 2
  • 10
  • 20