Comment vectoriser une simulation de marche aléatoire dans Matlab

Question

Je suis en train de réécrire un modèle de simulation de monte carlo dans Matlab en mettant l'accent sur la lisibilité. Le modèle implique de nombreuses particules (représentées par (x, y, z)) suite à une marche aléatoire sur un petit ensemble d'états avec certaines probabilités de terminaison. L'information pertinente pour la sortie est le nombre de particules qui se terminent dans un état donné.

La simulation nécessite suffisamment de particules pour que chaque particule soit individuellement coûteuse. La vectorisation semble être le moyen d'obtenir des performances de Matlab, mais existe-t-il une manière idiomatique de créer une version vectorisée de cette simulation dans Matlab? Je me bats la tête contre le mur pour y arriver - j'ai même essayé de créer une matrice nStates x nParticles représentant chaque combinaison particule-état, mais cette approche devient rapidement hors de contrôle (lisibilité) car les particules rebondissent état à déclarer indépendamment l'un de l'autre. Devrais-je juste mordre la balle et passer à une langue plus appropriée pour cela?

Answer 1

Il suffit d'écrire le code comme vous le feriez normalement. Presque toutes les fonctions de matlab peuvent accepter et retourner une entrée vectorisée. Par exemple, pour simuler un mouvement brownien de particules N en 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 

si vous voulez avoir un autre sigma pour chaque position, vous feriez sigma un vecteur de la même taille que la position et l'utilisation « des temps de point » notation pour indiquer l'élément par élément opération

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

si la diffusion était une fonction arbitraire de position et un élément aléatoire, vous devriez juste écrire une fonction vectorisée, par exemple

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

et ainsi de suite