Trackpy与手动计算布朗粒子MSD的差异原因排查

我正在研究trackpy算法与手动方法计算均方位移（MSD）的差异，计划用trackpy分析视频中布朗扩散粒子的轨迹。

1. 生成2D布朗扩散粒子轨迹 #

使用以下代码生成轨迹：

# time interval for steps (comparable to camera fps)
dt = 1/485                             
# total measurement time (comparable to video length)
totalT = 4850                         
# Theoretical diffusion constant for a particle of size 30nm
diffC = 16e-12           
              
track = np.zeros((totalT,3))

for i in range(1,len(track)):
    track[i,0] = track[i-1,0]+dt
    track[i,1] = track[i-1,1] + np.sqrt(2*diffC*dt)*np.random.normal(0,1)
    track[i,2] = track[i-1,2] + np.sqrt(2*diffC*dt)*np.random.normal(0,1) 

2. 手动计算MSD #

接着用以下代码手动计算所有时间间隔的MSD：

MSD = np.zeros((len(track),4))

for tau in range(1,len(track)):
    step = 1
    displacements = track[0:-tau:step,1:]-track[tau::step,1:]
    MSD[tau,0] = tau*dt
    MSD[tau,1] = np.mean(displacements[:,0]**2)-np.mean(displacements[:,0])**2
    MSD[tau,2] = np.mean(displacements[:,1]**2)-np.mean(displacements[:,1])**2
    MSD[tau,3] = MSD[tau, 1] + MSD[tau, 2]

3. trackpy计算MSD #

之后使用tp.motion.msd通过trackpy计算MSD。

实验结果 #

• 调整trackpy的max_lagtime参数后，MSD与时间延迟呈线性关系，得到的扩散常数与生成轨迹时的设定值一致。
• 手动计算的MSD与时间延迟无线性关系，无法得到扩散常数。

我的疑问：两种方法使用的数据点数量完全相同，为何手动计算无法得到线性关系？是trackpy算法存在问题，还是我的手动计算逻辑有误？

注：已查阅trackpy文档，了解到无丢帧情况下它使用FFT计算MSD，这是更快的方法，但我认为不会影响准确性。

你的手动计算逻辑存在错误，问题出在MSD的计算公式上。

对于布朗运动的MSD，正确定义是位移平方的平均值，即E[(x(t+τ) - x(t))²]，但你手动计算时额外减去了(np.mean(displacements))²——这是计算方差的操作，完全不符合MSD的定义。

错误原因 #

布朗运动中粒子的位移Δx = x(t+τ) - x(t)本身是无定向的，其平均值理论上为0，此时E[Δx²]等于位移的方差。但当τ增大时，可用的位移样本数减少，np.mean(displacements)的随机波动会被放大，你额外减去的平方项会干扰MSD的计算结果，导致曲线偏离线性。

修正后的手动计算代码 #

MSD = np.zeros((len(track),4))

for tau in range(1,len(track)):
    step = 1
    displacements = track[tau::step,1:] - track[0:-tau:step,1:]
    MSD[tau,0] = tau*dt
    MSD[tau,1] = np.mean(displacements[:,0]**2)  # 直接计算位移平方的均值
    MSD[tau,2] = np.mean(displacements[:,1]**2)
    MSD[tau,3] = MSD[tau, 1] + MSD[tau, 2]

补充说明 #

trackpy的tp.motion.msd严格遵循MSD的定义计算，FFT方法只是加速计算的手段，结果和直接求和的精确计算完全一致，因此能得到符合预期的线性关系。

内容的提问来源于stack exchange，提问作者nithilamk