Расчет среднеквадратичного смещения (MSD) с MATLAB

Question

Это часть вступительного курса Matlab для биологов.У меня есть точки данных (для одной частицы!) В матрице с 4 столбцами (время, x, y, z) и несколькими тысячами строк.Я хочу вычислить среднеквадратичное смещение для частицы, используя координаты XYZ для всех временных шагов.MSD определяется как MSD = среднее (r (t) -r (0)) ^ 2, где r (t) - положение частицы в момент времени t, а r (0) - начальное положение, поэтому в некотором смысле эторасстояние, пройденное частицей за промежуток времени t.Это то, что я до сих пор.Я, очевидно, новичок в MATLAB, и я очень признателен за помощь / ввод сома.

dx=zeros(length(data),1);   %Create space for displacements.
dy=zeros(length(data),1);
dz=zeros(length(data),1);

for i=1:length(pos) 
    dx(i)=data(i+1,2)-data(1,2);   %Calculate the distance at each time step
    dy(i)=data(i+1,3)-data(1,3);   %back to the origin.
    dz(i)=data(i+1,4)-data(1,4);
end

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

примечание: dt = delta t= шаг по времени, Pos = позиция

• Рассчитать смещение между Pos 2 и Pos 1 (dt = 1)
• Рассчитать смещение между Pos 3 и Pos 1 (dt = 2)
• Рассчитать смещение между Pos 4 и Pos 1 (dt = 3) и т. Д. ...

Теперь вместо этого сравните каждую позицию с Pos 2.

• Рассчитать смещение между Pos 3 и Pos 2 (dt = 1)
• Вычислить смещение между Pos 4 и Pos 2 (dt = 2)
• Вычислить смещение междуPos 5 и Pos 2 (dt = 3) и т. Д. ...

Теперь сравните каждую позицию с Pos 3

• Рассчитайте смещение между Pos 4 и Pos 3 (dt = 1)
• Рассчитать смещение между Pos 5 и Pos 3 (dt = 2)
• Рассчитать смещение между Pos 6 и Pos 4 (dt = 3) и т. д ...

Положение 4 ...

• Рассчитать смещение между Pos 5 и Pos 4 (dt = 1)
• Рассчитать смещение между Pos 6 и Pos 4 (dt= 2) ... и т. Д. Для тысяч строк.

Таким образом, для каждой точки будет намного больше точек данных.Имеет ли смысл так поступать с кем-либо еще?Не сделает ли это вычисление менее «шумным» (т. Е. Сделает график msd vs более плавным?) Возможно, есть некоторые удобные функции matlab, чтобы облегчить эту проблему?Я бы очень признателен за какой-то вклад.

Answer 1

Если вы вычисляете смещение только относительно первой позиции, то вы фактически ничего не усредняете, поскольку MSD (1) - это среднее расстояние, которое ваш объект проходит за один временной шаг.Таким образом, ваши рассуждения на месте.

Однако я бы рассчитал MSD наоборот: получить все смещения при dt = 1 (1-2,2-3,3-4, ...) и среднем.Это MSD (1).Затем вы получите все смещения при dt = 2 (1-3,2-4, ...) и в среднем.Это MSD (2).И так далее.

Полезным свойством Matlab является то, что вы можете векторизовать некоторые вычисления, т.е. выполнять вычисления для всего массива, а не делать их поэлементно.Поэтому, если у вас есть массив a с координатами 100 на 1, разница между каждой координатой и следующей будет b=a(2:100)-a(1:99), или, в более общем смысле, b=a(2:end)-a(1:end-1), поэтому b(1) равно a(2)-a(1), b(2) is a(3)-a(2) etc.

Чтобы вычислить MSD из вашего массива data (где я предполагаю, что время в равных шагах!), вы должны написать

nData = size(data,1); %# number of data points
numberOfdeltaT = floor(nData/4); %# for MSD, dt should be up to 1/4 of number of data points

msd = zeros(numberOfDeltaT,3); %# We'll store [mean, std, n]

%# calculate msd for all deltaT's

for dt = 1:numberOfDeltaT
   deltaCoords = data(1+dt:end,2:4) - data(1:end-dt,2:4);
   squaredDisplacement = sum(deltaCoords.^2,2); %# dx^2+dy^2+dz^2

   msd(dt,1) = mean(squaredDisplacement); %# average
   msd(dt,2) = std(squaredDisplacement); %# std
   msd(dt,3) = length(squaredDisplacement); %# n
end

Answer 2

Обратите внимание, что для анализа MSD доступен класс MATLAB: MSD анализатор , включая обширную документацию и учебные пособия .

Основное использование:

ma = msdanalyzer(2, SPACE_UNITS, TIME_UNITS);
ma = ma.addAll(tracks);
ma = ma.computeMSD;