# Help with speeding up my code and adding a 4 plot sub function/local function?

1 view (last 30 days)
Anon Anonson on 22 Feb 2014
So I have a code that is a Monte Carlo simulation of Buffon's needle experiment and I'm having some issues with it. I'm not terribly great at using Matlab, hence why I'm coming here for this. When I run the code, it takes hours to run but I know that there must be a way to speed it up, however, I have not found this way. Here is the code that I have thus far:
% This m-file (buffon.m) will produce a Monte Carlo simulation in which a
% distribution of estimates for the value of pi or trials will be found by
% completing a large number of separate needle drops or experiments known
% as Buffon's Needle Experiment. The number of experiments will then be
% increased until the average of the distribution of the number of trials
% will be within a given tolerance.
function [est,num_trials,num_exp] = buffon(num_trials,num_exp,tol)
est = 0;
while abs(pi-est) > tol
if num_exp > 10000
disp('Max allowable experiments exceeded')
break
end
for i = 1:num_trials
center_mass = rand(1,2,num_exp);
angle = rand(1,num_exp)*pi;
cross = 0;
for j = 1:num_exp
if center_mass(1,2,j) + 0.5*sin(angle(j)) > 1 || center_mass(1,2,j) - 0.5*sin(angle(j)) < 0
cross = cross + 1;
end
end
prob = cross/num_exp;
estimate_pi(i) = 2/prob;
end
est = mean(estimate_pi);
num_trials = num_trials*10;
end
end
1. As for the issue of how long this takes, I have no ideas on how to speed this code up. But what I think I've found is that the part that takes the longest is this one:
for i = 1:num_trials
center_mass = rand(1,2,num_exp);
angle = rand(1,num_exp)*pi;
and I believe it is the
rand(num_exp)*pi
that is taking forever.
Does anyone know of any ways to possibly speed this up at all? Any advice would help.
2. My next question is this: I want to add in a part at the end of the code to include a separate local function that plots four separate visual figures that include the following:
-Figure 1: a scatter plot representation of the ability of the function rand to provide an unbiased input using the data associated with the angle of the needle. This plot should show the angle as a function of the experiment for a given trial. -Figure 2: a cartoon representation (MATLAB-generated plot) of the distribution of needles for one trial plotted on a square containing marked vertical separations are distinguishable from those needles that do not. This plot should represent the needles as lines of different colors depending on whether they do or do not cross a vertical separating line, and the plot should also show the vertical separating lines on the square playing field. -Figure 3: a histogram (using the function histfit) binning the estimated values pi for each trial (or for each set of experiments), quoting the mean and the standard deviation of the set of trials in the legend or in a text label. -Figure 4: a scatter plot showing how the common logarithm (i.e., log10) of the absolute value of the error in the estimate varies as a function of the common logarithm of the number of experiments completed.
Sammy Williamson on 9 Aug 2019
function approxPi = buffons(totalNeedles)
numberCrossing = zeros(1, totalNeedles+1);
for i = 1:totalNeedles
xOffset = rand()*0.5;
angle = rand() * pi / 2;
height = 0.5 *sin(angle);
if xOffset <= height
numberCrossing(i+1) = numberCrossing(i) + 1;
else
numberCrossing(i+1) = numberCrossing(i);
end
end
approxPi = 2*(1:totalNeedles)./numberCrossing(2:end);
end

Nitin on 22 Feb 2014
Have you considered using parfor instead of a for loop?