matlab - 蒙特卡罗方法的向量化来近似 pi

Question

我成功地编写了以下代码，基于 for 循环，使用蒙特卡罗方法近似数字 pi:

function piapprox = calcPiMC(n)
count = 0;      % count variable, start value set to zero
for i=1:n;      % for loop initialized at one
    x=rand;     % rand variable within [0,1]
    y=rand;     

    if (x^2+y^2 <=1)    % if clause, determines if pair is within the first quadrant of the unit circle
    count = count +1;   % if yes, it increases the count by one
    end

end

piapprox = 4*(count/n);

end

这个方法很棒，但是对于较大的 n 值来说会很困难。作为我自己的一项任务，我想尝试对以下代码进行矢量化，不使用任何循环等。在我的脑海中，我可以想到一些对我来说有意义的伪代码，但到目前为止我还无法完成它。这是我的方法:

function piapprox = calcPiMCVec(n)

count = zeros(size(n));    % creating a count vector filled with zeros
x = rand(size(n));         % a random vector of length n
y = rand(size(n));         % another random vector of length n

if (x.*x + y.*y <=  ones(size(n)))   % element wise multiplication (!!!)
    count = count+1;                 % definitely wrong but clueless here
end

piapprox=4*(count./n); 
end 

我希望在阅读这个伪代码时我的想法看起来很清楚，但是我会对它们进行评论。

• 我开始创建一个由零组成的计数向量，其长度由向量 n 的大小决定
• 接下来，我对随机条目 x 和 y 做了完全相同的事情
• 与填充相同长度的向量相比，if 子句应该确定哪些条目小于 1，但我非常确定这段代码是错误的。
• 最后我想将 if 子句为真的向量的数量存储到一个向量中，这段代码也是错误的，我想我必须在这里使用 cumsum 函数，但结果是错误的

Answer 1

这是一个矢量化解决方案，其中有一堆评论，部分引用了您的尝试。

function piapprox = calcPiM(n)
%#CALCPIM calculates pi by Monte-Carlo simulation
%# start a function with at least a H1-line, even better: explain fcn, input, and output

%# start with some input checking
if nargin < 1 || isempty(n) || ~isscalar(n)
   error('please supply a scalar n to calcPiM')
end

%# create n pairs of x/y coordinates, i.e. a n-by-2 array
%# rand(size(n)) is a single random variable, since size(n) is [1 1]
xy = rand(n,2); 

%# first, do element-wise squaring of array, then sum x and y (sum(xy.^2,2))
%# second, count all the radii smaller/equal 1
count = sum( sum( xy.^2, 2) <= 1);

%# no need for array-divide here, since both count and n are scalar
piapprox = 4*count/n;