Cómo vectorizar una simulación de paseo aleatorio en Matlab

Question

estoy reescribiendo un modelo de simulación Monte Carlo en Matlab con un énfasis en la legibilidad. El modelo implica muchas partículas (representados como (x, y, z)) después de un paseo aleatorio sobre un pequeño conjunto de estados con ciertas probabilidades de terminación. La información relevante para la salida es el número de partículas que terminan en un estado dado.

La simulación requiere suficientes partículas que se ejecuta para cada partícula es individualmente prohibitivo coste. Vectorización parece ser la forma de obtener rendimiento de Matlab, pero ¿hay alguna manera idiomática de crear una versión vectorizada de esta simulación en Matlab?

Estoy batiendo la cabeza contra la pared para lograr esto - Incluso he intentado crear un nStates x nParticles matriz que representa cada combinación de estado de partículas, pero este enfoque rápidamente se sale de control (lectura) ya que las partículas rebotan de un estado afirmar independientemente uno de otro. ¿Debo de tripas corazón y cambio a un lenguaje más adecuado para esto?

Answer 1

Sólo tiene que escribir el código como lo haría normalmente. Casi todas las funciones de MATLAB pueden aceptar y devolver la entrada vectorizado. Por ejemplo, para simular un movimiento browniano de partículas N en 1 dimensión

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

si quería tener un sigma diferente para cada posición, que haría sigma un vector del mismo tamaño que la posición y el uso "tiempos punto" notación para indicar elemento por elemento de operación

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

Si la dispersión es una función arbitraria de posición y algún elemento aleatorio, sólo tendría que escribir una función vectorizado, por ejemplo.

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);

y así sucesivamente