Optimizing Monte Carlo Simulation

I am running into issues with the simulation of 1 million particles or more. I need the code to run under a minute at 10 million particle sims. I estimate it runs at about 22min for 10 million particles currently. If anyone could help me resolve this or get the computational time down, it would be great. I have the parallel computing toolbox downloaded if it is neccesary.

clear;clc 
tic
N = 10000000;
EndPos = zeros(3,N);
SigmaA = 0.00032;
SigmaS = .35429;
Sigmat = SigmaA + SigmaS;
Pscat = SigmaS / Sigmat;
Pabs = SigmaA / Sigmat;
% Precompute random angles
Thetas = unifrnd(0, 2*pi, 1, N);
Phis = unifrnd(0, pi, 1, N);
DistProbs = unifrnd(0, 1, 1, N);
Distance = -1 / Sigmat * log(1 - DistProbs);
XYZ = [Distance .* sin(Thetas) .* cos(Phis); Distance .* sin(Thetas) .* sin(Phis); Distance .* cos(Thetas)];
InteractionProbs = unifrnd(0, 1, 1, N);
ScatEv = InteractionProbs > Pabs;
ScatDist = XYZ(:,ScatEv);
AbsDist = XYZ(:,~ScatEv);
EndPos(:,1:width(AbsDist)) = AbsDist;
while ~isempty(ScatEv)
    Thetas = unifrnd(0, 2 * pi, 1, width(ScatDist));
    Phis = unifrnd(0, pi, 1, width(ScatDist));
    DistProbs = unifrnd(0, 1, 1, width(ScatDist));
    Distance = -1 / Sigmat * log(1 - DistProbs);
    XYZi = [Distance .* sin(Thetas) .* cos(Phis); Distance .* sin(Thetas) .* sin(Phis); Distance .* cos(Thetas)];
    InteractionProbs = unifrnd(0, 1, 1, width(ScatDist)); 
    ScatEv = InteractionProbs > Pabs;
    ScatDisti = XYZ(:,ScatEv);
    AbsDisti = XYZ(:,~ScatEv);
    ScatDisti = ScatDist(:, ScatEv) + ScatDisti;
    AbsDist = ScatDist(:, ~ScatEv) + AbsDisti;
    ScatDist = ScatDisti;
    EndPosInd = find(EndPos(1,:) == 0, 1,'first');
    EndPosInd1 = EndPosInd + width(AbsDist) - 1;
    EndPos(:,EndPosInd:EndPosInd1) = AbsDist;
    
end
PosSquare = EndPos.^2;
rsquare = sum(PosSquare);
rsquarebar = sum(rsquare)/N;
RquareSTD = std(rsquare);
RsquareError = RquareSTD/sqrt(width(rsquare));
Lsquare = 1/6*rsquarebar; %diffusion area estimate
toc

Answer 1

[Walter Roberson]

You can make some micro-speedups.

A = unifrnd(0, 2*pi, 1, N)

is slightly slower than

A = 2 * pi * rand(1, N);

width(X) is slightly slower than size(X,2)

Distance = -1 / Sigmat * log(1 - DistProbs);

would be slightly faster as

RNSigmat = -1 ./ Sigmat;

before the loop and then

Distance = log(1 - DistProbs) .* RNSigmat;

but I would not expect those changes to add up to even 1 second savings in execution time.

[Sam Marshalik]

I am not sure if Parallel Computing Toolbox will be helpful with this particular problem, but you could run multiple simulations at the same time. Meaning, I do not think Parallel Computing Toolbox will help speed up a single simulation but you could run multiple simulations in parallel at the same time on your machine with parfor.