Higher-dimensional Monte Carlo Integration

HigherMonteCarlo.py

#
# author Huber Flores
#

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)