Comment vectoriser une simulation de marche aléatoire dans Matlab

Question

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

La simulation nécessite suffisamment de particules pour que l'exécution de chaque particule est individuellement prohibitif des coûts. Vectorisation semble être le moyen d'obtenir des performances de Matlab, mais est-il un moyen idiomatiques de créer une version vectorisée de cette simulation en Matlab?

Je battre ma tête contre le mur pour y arriver - j'ai même essayé créer un nStates x matrice nParticles représentant chaque combinaison de particules d'état, mais cette approche des spirales rapidement hors de contrôle (lisibilité) depuis rebond des particules d'un état d'indiquer indépendamment l'un de l'autre. Dois-je mordre juste la balle et passer à un langage plus approprié pour cela?

Answer 1

Il suffit d'écrire le code comme vous le feriez normalement. Presque toutes les fonctions de Matlab peuvent accepter et renvoyer entrée vectorisé. Par exemple, pour simuler un mouvement brownien des particules dans une dimension 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 vous voulez avoir un autre sigma pour chaque position, vous feriez sigma un vecteur de la même taille que la position et utiliser la notation « dot fois » pour indiquer l'élément par élément de fonctionnement

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