I would like to know if I calculate the MSD correctly, here is a Jupyter snippet I am using:

import numpy as np
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
%matplotlib inline

x = np.random.randn(100)

df = pd.DataFrame(x, columns=['x'])
df['time'] = np.arange(0,len(df))
df['MSD'] = pd.Series((x - x[0])**2).rolling(window=10, min_periods=1).mean()

x_plt, = plt.plot(x, label='x(t)')
msd_plt, = plt.plot(df['MSD'], label='MSD(t)')
plt.legend(handles=[line_up, line_down])

enter image description here

Is there a smart way to average over the variances? The original definition of MSD has average inside, but I don't relay averaging here because I have only one particle.


1 Answer 1


According to the article on which you base your MSD computation, it's specified as a single number (an average over N items) which is defined as the mean (average) of the sum of squared-distances between positions of $N$ particles at time $t$ and given reference positions $x_n(0)$ for each particle at a reference time.

Assuming that you want to compute MSD over a 1D particle displacement of $N$ items, you would use the following Matlab/Octave code to achieve it:

N = 1024;                 % Number of particles
std1 = 2.34;              % set some variance
xpos = std1 *sqrt(3)*2*( randn(1,N) - 0.5); % N particles with positions. 
xref = 0;                 % using same reference position for each particle.
MSD = (1/N)* sum ( (xpos-xref).^2 );  % compute MSD of N particles at time t

figure,plot(xpos);        % plot positions of each particle at time t
hold on;                  % plot MSD of those N particles at time t.
plot(MSD*ones(1,N));      % MSD*ones(1,N) is used to plot N points for MSD
title('MSD over particle positions at time t');

Assuming that you also want to plot MSD as a function of sampled time, you could then then use the following code to achieve it:

N = 1024;                 % Number of particles
M = 128;                  % Number of MSD computations sample time points
std1 = 2.34;              % Set some variance, this can change over time.
xpos = std1 *sqrt(3)*2*( randn(N,M) - 0.5); % N particles for each M time samples. 
xref = 0;                 % using same reference position for all time and each poarticle.
MSD = (1/N)* sum ( (xpos-xref).^2 );  % now MSD is an 1 x M vector

figure,plot(MSD);        % plot particle positions and MSD
  • $\begingroup$ Why does N=1024? I am not sure how you adjusted the algorithm to one particle. $\endgroup$
    – 0x90
    Jul 12, 2017 at 10:48
  • $\begingroup$ N is the number of particles in the computation, you can set it according to your problem. $\endgroup$
    – Fat32
    Jul 12, 2017 at 10:51
  • $\begingroup$ so for single particle N=1. $\endgroup$
    – 0x90
    Jul 12, 2017 at 10:54
  • 1
    $\begingroup$ yes for a single particle N=1. $\endgroup$
    – Fat32
    Jul 12, 2017 at 11:00

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.