用MATLAB计算均方位移(msd)

Question

这是生物学家matlab入门课程的一部分.我有一个矩阵的数据点(对于单个粒子!)有4列(时间,x,y,z)和几千行.我想要做的是使用所有时间步长的xyz坐标计算粒子的均方位移.MSD定义为MSD = average(r(t)-r(0))^ 2其中r(t)是粒子在时间t的位置,r(0)是初始位置,所以在某种意义上它是粒子在时间间隔t内行进的距离.这就是我到目前为止所拥有的.我显然是一个MATLAB新手,我真的很感激som帮助/输入.

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)...等数千行.

这样,每个dt会有更多的数据点.对其他人这样做是否有意义？这会使计算减少"嘈杂"(即使msd与时间的关系曲线更平滑吗？)也许有一些方便的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的一个有用属性是你可以对一些计算进行矢量化,即对整个数组进行计算,而不是逐个元素地进行计算.因此,如果你有一个数组a100乘1坐标,坐标和下一个是各自之间的差值b=a(2:100)-a(1:99),或者更一般地,b=a(2:end)-a(1:end-1),所以b(1)就是a(2)-a(1),b(2)是a(3)-a(2)等

要从你的数组计算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