在 Julia 中使用 Monte-Carlo 方法计算 n-ball 的体积
Calculating volume of n-ball using Monte-Carlo method in Julia
我想写一个计算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
我不确定出了什么问题
好的,所以有一些突出的错误。您的 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
和一个函数,这样它就不必为方块分配中间数组。我们还可以使用 @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.
我想写一个计算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
我不确定出了什么问题
好的,所以有一些突出的错误。您的 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
和一个函数,这样它就不必为方块分配中间数组。我们还可以使用 @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.