将蒙特卡罗方法向量化为近似 pi

Question

我已经成功编写了以下代码，基于 for 循环使用 Monte-Carlo 方法逼近数字 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

这是一个矢量化解决方案，其中包含一些部分涉及您的尝试的 cmets。

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;