I want to use the Monte-Carlo integration method to find the volume of an n-sphere and compare it to the analytical result. I want to do this for 10^6 points and 10^9 points and while it works for 10^6 points somewhat (takes about a minute for n = 2 (circle) , n = 3 (sphere) and n = 12), it is awfully slow for 10^9 points.
Short explanation of MC-Method: To find the Volume of a n-Sphere with radius r = 1, I imagine an easy known Volume (say a n-cube with sides of length 2*r) which fully contains the n-sphere. Then I sample from a uniform distribution points in the n-cube and check, if the point lies in the sphere. I count all the so generated points inside the n-sphere. The ratio of V_sphere/V_cube can be approximated as N_inside/N_total and therefore, V_sphere = V_cube * N_inside/N_total
Here is the function:
def hyp_sphere_mc(d,samples):
inside = 0 #number of points inside sphere
sum = 0 #sum of squared components
for j in range(0,samples):
x2 = np.random.uniform(0,1,d)**2 #generate random point in d-dimensions
sum = np.sum(x2) #sum its components
if(sum < 1.0):
inside += 1 #count points inside sphere
V = ((2)**d)*(float(inside)/samples) #V = V_cube * N_inside/N_total
V_true = float(math.pi**(float(d)/2))/math.gamma(float(d)/2 + 1) #analytical result
ERR = (float(abs(V_true-V))/V_true)*100 #relative Error
print "Numerical:", V, "\t" , "Exact: ", V_true, "\t", "Error: ", ERR
I guess the problem is, that for every iteration, I generate a new random array and this takes a lot of time, especially if I have 10^9 iterations. Is there any way to speed this up?
You can replace the loop with the following:
inside = np.sum(np.sum(np.random.rand(samples,d)**2,1)<1)
When using numpy you should try to avoid loops. The idea is that you can just generate all the samples at once in a matrix, and then vectorize all subsequent operations.