julia - 在 Julia 中使用 Monte-Carlo 方法计算 n-ball 的体积

Question

我想编写一个代码来计算 n 球的体积。
使用一系列随机点 r_i(i = 1, ... , N) 可以估算出这种球的体积:
其中d是维度，N是点数

我在运行这段代码时遇到问题

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

我不知道哪里出了问题

Answer 1

好的，所以有一些突出的错误。您的 end 比您需要的多了一个，并且行

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

有几个明显的错误，因为 (1) 您在 julia 中使用 [] 索引数组，而不是 ()，以及 (2) rand(n) 生成 n 个随机数，而不是 0 和 n 之间的随机数。所以如果我们解决这些问题

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

现在我们有了一种方法，如果我们测试简单的 3d 情况，它似乎可以为我们提供正确的数字

julia> MonteCarloHypersphereVolume(1,3,100000)
4.18328

julia> 4/3*pi*1
4.1887902047863905

现在，我们可以进行一些优化。目前这需要大约 13 毫秒和大量分配

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.

我们可以做得更好，方法是将 randoms = zeros( 1, dimension ) 移到循环之外，并将 sum 与一个函数一起使用，这样它就不必为方 block 分配一个中间数组。我们还可以使用 @inbounds 来加快 for 循环的速度，因为我们不应该期望这里有任何越界索引。

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

我也做了一些风格上的改变。除其他事项外，在 Julia 中，函数通常按照约定使用小写字母(而类型可以是 CamelCase)。

无论如何，这给了我们大约 6 倍的加速，并显着减少了内存使用和分配

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.