Higher-dimesional Monte Carlo Integration (Python)

#parameters
* fd(x,y,z) is x*x+y*y+z*z<=1
* f(x,y,z) is 1.0
* box is [(-1.0,1.0), (-1.0,1.0), (-1.0,1.0)]

should result in 4*Pi/3


import numpy as np
from scipy import *
import random

box = ([-1.0,1.0],[-1.0,1.0],[-1.0,1.0])

def volume(B):
vol = zeros(len(box))
for v in range(len(box)):
b = box[v]
vol[v] = b[1] - b[0]

volume = np.multiply.reduce(vol)
return volume

def fd(x,y,z):
if (x*x+y*y+z*z<=1):
return True
else:
return False

def f(x,y,z):
return 1.0

#generating values
dim = 3
points = 1000
values = zeros(dim*points)
for i in range(dim*points):
values[i] = random.uniform(-1,1)

#xyz values vector
vector = values.reshape(points,dim)

#vectorizing functions
vfd = np.vectorize(fd)
vf = np.vectorize(f)

#calculating values of the vector using fd
fd_values = zeros(points)
for j in range(points):
fd_values[j] = vfd(vector[j,0],vector[j,1],vector[j,2])

sum = 0
for h in range(points):
if (fd_values[h]==1):
sum = sum + vf(vector[h,0],vector[h,1],vector[h,2])

#Alternative
#print fd_values
#print sum(where(fd_values>0,f(x,y,z),0))/float(dim)

approx = (volume(box)/float(points))*sum

print "Multidimensional Integration is:" + str( approx)



This is a working document. If there are any errors, or if you disagree with anything which I have said, or have any suggestions for improvements please email me kyoraul@gmail.com.