[解決済み] MATLABによる平均二乗変位(msd)の算出

質問内容

これは生物学者のためのmatlab入門コースの一部です。私は4列（時間、x、y、z）と数千行の行列にデータポイント（単一の粒子！）を持っています。私がしたいことは、すべてのタイムステップのxyz座標を使用して、粒子の平均二乗変位を計算することです。MSDはMSD=average(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 = デルタ t = 時間ステップ、Pos = 位置

• Pos2とPos1の間の変位を計算する(dt=1)
• Pos 3とPos 1の間の変位を計算する(dt=2)
• Pos4とPos1間の変位(dt=3)を計算する...など。

ここで、代わりにすべての位置をPos 2と比較します。

• Pos 3とPos 2の間の変位を計算する(dt=1)
• Pos4とPos2間の変位量を計算する(dt=2)
• Pos 5とPos 2の間の変位を計算する(dt=3)...といった具合に。

ここで、各ポジションをPos 3と比較します

• Pos 4 と Pos 3 の間の変位を計算する (dt=1)
• Pos5とPos3間の変位量を計算する(dt=2)
• Pos6とPos4間の変位量(dt=3)を計算する...といった具合です。

ポジション4...

• Pos5とPos4の間の変位を計算する(dt=1)
• Pos 6 と Pos 4 の間の変位を計算 (dt=2) ...等、数千行分。

この方法では、各dtに対してより多くのデータポイントが存在することになります。このようにすることに意味があるのでしょうか？このようにすれば、計算の「ノイズ」が少なくなるでしょうか（つまり、msd対時間のプロットが滑らかになるでしょうか）。もしかしたら、この問題を簡単にする便利なmatlab関数があるかもしれませんね。私は本当にいくつかの意見に感謝します。

どのように解決するのですか？

最初の位置に対してのみ変位を計算するのであれば、MSD(1)は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) は 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