Come per vettorizzare una simulazione random walk in Matlab

Question

Sto riscrivendo un monte carlo modello di simulazione in Matlab con l'accento sulla leggibilità. Il modello prevede molte particelle (rappresentati come (x, y, z)) seguendo un cammino casuale su un piccolo insieme di stati determinati probabilità di terminazione. Le informazioni rilevanti per l'uscita è il numero di particelle che terminano in un dato stato.

La simulazione richiede abbastanza particelle che esegue per ogni particella è singolarmente costo proibitivo. Vettorializzazione sembra essere il modo per ottenere le prestazioni di Matlab, ma c'è un modo idiomatico di creare una versione di vettorializzare questa simulazione in Matlab?

sto battendo la testa contro il muro per raggiungere questo obiettivo - ho anche cercato di creare un nStates x nParticles matrice che rappresenta ogni combinazione particella-stato, ma questo approccio spirali rapidamente fuori controllo (leggibilità) poiché le particelle rimbalzano da stato affermare indipendentemente l'uno dall'altro. Devo solo stringere i denti e passare a un linguaggio più adatto a questo?

Answer 1

Basta scrivere il codice come si farebbe normalmente. Quasi tutte le funzioni MATLAB può accettare e tornare in ingresso vettorializzare. Ad esempio, per simulare un moto browniano delle particelle N in 1 dimensione

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

se si voleva avere un sigma diverso per ogni posizione, si farebbe sigma un vettore della stessa dimensione come la posizione e l'uso "dot volte" notazione per indicare elemento per effetto elemento

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

se la dispersione è una funzione arbitraria della posizione e qualche elemento casuale, sarebbe solo devono scrivere una funzione vectorized, per esempio.

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

e così via