Как векторизовать симуляцию случайного блуждания в Matlab - PullRequest
Новая
       48

Как векторизовать симуляцию случайного блуждания в Matlab

3 голосов
/ 24 сентября 2010

Я переписываю имитационную модель Монте-Карло в Matlab с акцентом на удобочитаемость. Модель включает в себя множество частиц (представленных в виде (x, y, z)), следующих за случайным блужданием по небольшому набору состояний с определенными вероятностями завершения. Информация, релевантная для вывода, представляет собой количество частиц, которые оканчиваются в данном состоянии.

Для моделирования требуется достаточное количество частиц, чтобы запустить его для каждой частицы в отдельности из-за высокой стоимости. Кажется, что векторизация - это способ повысить производительность Matlab, но есть ли какой-то идиоматический способ создания векторизованной версии этого моделирования в Matlab?

Я бьюсь головой об стену, чтобы добиться этого - я даже пытался создать матрицу nStates x nParticles, представляющую каждую комбинацию состояния частицы, но этот подход быстро выходит из-под контроля (читаемость), так как частицы отскакивают от состояния заявлять независимо друг от друга. Стоит ли просто прикусить пулю и переключиться на язык, более подходящий для этого?

1 Ответ

3 голосов
/ 24 сентября 2010

Просто напишите код, как обычно. Почти все функции Matlab могут принимать и возвращать векторизованный ввод. Например, для моделирования броуновского движения N частиц в 1 измерении 

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

если бы вы хотели иметь разные сигмы для каждой позиции, вы бы сделали сигма вектором того же размера, что и положение, и использовали бы нотацию "точек раз", чтобы указать элемент при операции элемента 

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

если бы рассеяние было произвольной функцией положения и некоторого случайного элемента, вам просто нужно написать векторизованную функцию, например 

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

и т. Д.

...