CODEWITHSUNDEEP
matlab
Spaced
Spaced

Reputation: 231

Vectorization of the Monte-Carlo-method to approximate pi

With success I have written the following code, based on a for loop to approximate the number pi using the Monte-Carlo-method:

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

This method is great but begins to struggle a lot for large values of n. As a task for myself I wanted to try to vectorize the following code, without using any loops or the like. In my head I can think of some pseudocode that would make sense to me, but I couldn't finish it so far. Here is my approach:

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

I hope my ideas seem clear when reading this pseudocode, I will however comment on them.

Upvotes: 2

Views: 1589

Answers (1)

Jonas
Jonas

Reputation: 74940

Here's a vectorized solution with a bunch of comments that partly refer to your attempt.

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;

Upvotes: 4

Related Questions