Higher-dimesional Monte Carlo Integration (Python)

	#
	# 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)

view raw HigherMonteCarlo.py hosted with ❤ by GitHub

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.