How to vectorize a random walk simulation in MATLAB

Question

I am rewriting a Monte Carlo simulation model in MATLAB with an emphasis on readability. The model involves many particles, represented as (x,y,z), following a random walk over a small set of states with certain termination probabilities. The information relevant for output is the number of particles that terminate in a given state.

The simulation requires enough particles that running it for each particle individually is cost prohibitive. Vectorization seems to be the way to get performance out of MATLAB, but is there any idiomatic way of creating a vectorized version of this simulation in MATLAB?

I'm beating my head against the wall to accomplish this - I've even tried creating a (nStates x nParticles) matrix representing each particle-state combination, but this approach quickly spirals out of control in terms of readability since particles bounce from state to state independently of one another. Should I just bite the bullet and switch to a language more suitable for this?

Witcher

Hard to make useful comments when you give so few details about the nature of the problem. – Nightrider

Apologies, I was trying to keep it abstracted for clarity. I am simulating light transport through layers of organic material. Each particle (photon) hits the material and can be reflected or refracted. This is determined with some probability proportional to material properties and angle of incidence. The particle can "bounce around" inside the layers. I want to know how many are reflected overall, and how many pass through. I do this by shooting many particles and tracing their paths. The question is, how to trace these paths as vectorized code rather than a while loop for each particle. – Witcher

Answer 1

Just write the code as you normally would. Almost all matlab functions can accept and return vectorized input. For instance, to simulate a brownian motion of N particles in 1 dimension

position = zeros([N 1]); %start at origin
sigma = sqrt(D * dt); %D is diffusion coefficient, dt is time step
for j = 1:numSteps
    position = position + sigma*randn(size(position));
end

if you wanted to have a different sigma for each position, you would make sigma a vector the same size as position and use "dot times" notation to indicate element by element operation

position = position + sigma.*randn(size(position));

if the scattering was an arbitrary function of position and some random element, you would just have to write a vectorized function, e.g.

function newstep = step(position)
%diffusion in a overdamped harmonic potential
newstep = -dt*k*position + D*randn(size(position));

for j = 1:numsteps; position = position + step(position);

and so on

Peal

This is really helpful. I have one additional question: in my case the particle is represented by velocity in 3D and the layer number (integer), i.e. a 4-tuple. I need different behavior depending on the layer number. How could we change the example to have multiple functions depending on the first element of the tuple and keep it vectorized? E.g. there could be 6 "newstep" functions, the choice of which depends on which layer. – Witcher

Can I use above function to simulate a set of mobile users. If yes how can include user parameters such as transmit power and interference into above function. – Discovert