Среднеквадратичное смещение (MATLAB)

Question

Пожалуйста, вы можете помочь мне понять, как рассчитать среднеквадратичное смещение для отдельной частицы, движущейся случайным образом в течение заданного периода времени. Я прочитал много статей по этому вопросу (включая Saxton, 1991, «Отслеживание отдельных частиц: распределение коэффициентов диффузии»), но все еще в замешательстве (не получил правильного ответа).

Позвольте мне начать с показа вам, как я это делаю, и, пожалуйста, поправьте меня, если я ошибаюсь: Я делаю это следующим образом:
1. Запустите программу с t = 0 до t = 100
2. Рассчитайте смещение (s (t) -s (t + tau)) на каждом временном шаге (т. Е. При t = 1,2,3, ... 100) и сохраните его в векторе
3. Квадратный ответ на номер 2
4. найти среднее к ответу 3

По сути, это то, что я делаю в Matlab

% Инициализируйте решетку квадратом, состоящим из 16 ненулевых узлов решетки, затем выполните% следующим образом для вычисления MSD:

for t=1:tend    
% Allow the particle to move randomly in the lattice. Then do the following    
[row,col]=find(lattice>0);    
centroid=mean([row col]);    
xvec=[xvec centroid(2)];    
yvec=[yvec centroid(1)];    

k=length(xvec)-1;  % Time    
dt=1;       
diffx = xvec(1:k) - xvec((1+dt):(k+dt));    
diffy = yvec(1:k) - yvec((1+dt):(k+dt));    

xsquare = diffx.^2;    
ysquare = diffy.^2;    
MSD=mean(xsquare+ysquare);    
end

Я пытаюсь найти MSD для вычисления коэффициента диффузии. Обратите внимание, что я моделирую коллекцию узлов решетки (16) для представления одной частицы (более биологически реалистичной) вместо одной. Я был кратким с комментарием в цикле for, так как он довольно длинный, но я рад отправить его вам.

Пока что я получаю очень маленькие значения MSD (в диапазоне 0,001-1), тогда как я должен получать значения в диапазоне (10-50). Частица движется на очень большие расстояния, поэтому мой диапазон 0,001-1 не может быть правильным!

---

Это выдержка из статьи, в которой я пытаюсь воспроизвести их фигуру:

"Мы начали с запуска некоторых симуляций в 1D для одного клетка. Мы позволили клетке двигаться на заданное количество Шаг по времени Монте-Карло (MCS), вычисленное среднее значение расстояние, пройденное за это время, повторило этот процесс 500 раз, и оцените среднеквадратичное расстояние для этого t. Затем мы повторили этот процесс десять раз, чтобы получить среднее значение , Причиной такого выбора повторений было сохранить время, необходимое для запуска моделирования в разумных пределах уровень, гарантирующий, что стандартное отклонение среднее было относительно мало (<7%) ". <br> Вы можете получить доступ к статье здесь "От дискретной к непрерывной модели биологического движения клеток, 2004, Turner et al., Physical Review E".

---

Любые советы приветствуются.

Answer 1

Прежде всего, почему частица покрывает несколько узлов решетки?Что в конечном итоге имеет значение для MSD, так это смещение центроида, которое можно представить в виде точки.Если ваша частица (или клетка) велика или имеет только большие шаги, вы всегда можете просто сделать более широкую сетку.Кроме того, если вы пытаетесь воспроизвести фигуру откуда-то еще, вам следует использовать тот же алгоритм.

Что вы делаете для симуляции Монте-Карло?Если все, что вам действительно нужно, это получить смещение, вы можете сгенерировать несколько случайных векторов движения за один раз (используя rand или randi) и использовать cumsum для вычисления позиций.Кроме того, вы планировали свои случайные прогулки, чтобы убедиться, что данные разумны?

Затем ваш код выглядит немного забавно (см. Комментарии).Почему бы вам просто не использовать код, указанный в этом ответе для расчета MSD по позициям?

for t=1:tend    
% Allow the particle to move randomly in the lattice. Then do the following    
[row,col]=find(lattice>0); %# what do you do this for?
centroid=mean([row col]);    
xvec=[xvec centroid(2)];    
yvec=[yvec centroid(1)];   %# till here, I have no idea what you want to do

k=length(xvec)-1;  % Time  %# you should subtract dt here
dt=1;                      %# dt should depend on t!
diffx = xvec(1:k) - xvec((1+dt):(k+dt));    
diffy = yvec(1:k) - yvec((1+dt):(k+dt));    

xsquare = diffx.^2;    
ysquare = diffy.^2;    
MSD=mean(xsquare+ysquare);    
end

Answer 2

На сколько размеров движется частица?

У меня сейчас нет Matlab, но вот как я могу это сделать в одном измерении:

% pos is the vector of positions
delta = pos(2:100) - pos(1:99);
meanSquared = mean(delta .* delta);