Calculating volume of n-ball using Monte-Carlo method in Julia

Question

I want to write a code that calculates the volume of a n-ball.
With a sequence of random points r_i(i = 1, ... , N) one can approximate the volume of such ball:
where d is the dimension and N the number of points

I am having trouble running this code

function MonteCarloHypersphereVolume(radius, dimension, number_of_generations)
    
number_within_sphere = 0;
  for i = 1 : number_of_generations
    randoms = zeros( 1, dimension );
    for j = 1 : dimension 
        randoms(j) = rand(radius * 2) - radius;
    end

    if sum( randoms .^ 2 ) <= radius^2
        number_within_sphere = number_within_sphere + 1;
    end
end

approximate_volume = (number_within_sphere/number_of_generations) * (2*radius)^dimension;

end

end

I am not sure what is going wrong

user_13639734 · Answer

Ok, so there are a few errors that stand out. You have one more end than you need, and the line

randoms(j) = rand(radius * 2) - radius;

has a couple apparent errors because (1) you index arrays with [] in julia, not (), and (2) rand(n) generates n random numbers, not random numbers between 0 and n. So if we fix those issues

function MonteCarloHypersphereVolume(radius, dimension, number_of_generations)
    
    number_within_sphere = 0
    for i = 1 : number_of_generations
        randoms = zeros( 1, dimension )
        for j = 1 : dimension 
            randoms[j] = 2*radius*rand() - radius
        end

        if sum( randoms .^ 2 ) <= radius^2
            number_within_sphere = number_within_sphere + 1
        end
    end

    approximate_volume = (number_within_sphere/number_of_generations) * (2*radius)^dimension;

end

Now we have an approach that seems to give us the right numbers if we test the simple 3d case

julia> MonteCarloHypersphereVolume(1,3,100000)
4.18328

julia> 4/3*pi*1
4.1887902047863905

Now, there are a few optimizations we can make. Currently this takes about 13 ms and a huge number of allocations

julia> using BenchmarkTools

julia> @benchmark MonteCarloHypersphereVolume(1,3,100000)
BechmarkTools.Trial: 375 samples with 1 evaluations.
 Range (min … max):  11.016 ms … 23.764 ms  ┊ GC (min … max): 0.00% … 12.36%
 Time  (median):     12.893 ms              ┊ GC (median):    0.00%
 Time  (mean ± σ):   13.339 ms ±  1.647 ms  ┊ GC (mean ± σ):  4.24% ±  4.51%

      ▄▃▂█  ▆▇▅▇▁▁ ▂▅ ▁                                        
  ▃▄█▄█████▇███████████▆▇▆▅▇▅█▅▅▄▃▆▃▄▄▄▃▁▃▃▃▁▁▃▄▃▁▃▁▁▁▃▁▁▃▁▃▃ ▄
  11 ms           Histogram: frequency by time        18.9 ms <

 Memory estimate: 21.36 MiB, allocs estimate: 200000.

We can do much better by moving randoms = zeros( 1, dimension ) outside the loop, and using sum with a function such that it doesn't have to allocate an intermediate array for the squares. We can also use @inbounds to speed up the for loops a bit given that we should not expect any out-of-bounds indices here.

function montecarlohyperspherevolume(radius, dimension, number_of_generations)
    
    number_within_sphere = 0
    randoms = zeros( 1, dimension )
    @inbounds for i = 1 : number_of_generations
        for j = 1 : dimension 
            randoms[j] = 2*radius*rand() - radius
        end

        if sum(abs2, randoms) <= radius^2 # could also use anonymous function x->x^2 if you didn't wnat to use abs2
            number_within_sphere += 1
        end
    end

    return approximate_volume = (number_within_sphere/number_of_generations) * (2*radius)^dimension

end

I've made a few stylistic changes as well. Among other things, in Julia functions are normally lower case by convention (while types can be CamelCase).

In any case, this gives us about a 6x speedup, and dramatically reduces memory usage and allocations

julia> @benchmark montecarlohyperspherevolume(1,3,100000)
BechmarkTools.Trial: 2249 samples with 1 evaluations.
 Range (min … max):  1.846 ms …   4.053 ms  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     2.080 ms               ┊ GC (median):    0.00%
 Time  (mean ± σ):   2.217 ms ± 377.123 μs  ┊ GC (mean ± σ):  0.00% ± 0.00%

   █▆▅▂▁▂                                                      
  ▅████████▇▆▅▅▅▃▃▃▃▃▃▄▃▃▂▃▂▃▂▂▃▂▂▂▂▂▁▁▂▁▂▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁ ▂
  1.85 ms         Histogram: frequency by time        3.61 ms <

 Memory estimate: 112 bytes, allocs estimate: 1.

Calculating volume of n-ball using Monte-Carlo method in Julia

1 Answers