1

I'd like to calculate the mean square displacement at all times of a 1d trajectory. I'm a novice programmer, so I've attempted it from a simulated random walk.

import numpy as np
import matplotlib.pyplot as plt

# generate random walk
path = np.zeros(60)
position = 0
for i in range(1, path.size):
    move = np.random.randint(0, 2)
    if move == 0: position += -1
    else: position += 1
    path[i] = position

# returns a vector of MSDs from a given path
def calcMSD(data):
    msd = np.zeros(data.size - 1)
    for i in range(1, data.size):
        squareddisplacements = (data[i:] - data[:data.size - i])**2
        msd[i - 1] = squareddisplacements.mean()
    return msd

msd = calcMSD(path)
plt.plot(msd)
plt.show()

Have I implemented this correctly? Any and all advice is appreciated, thanks.

Luke
  • 11
  • 2
  • what do you mean by "walk"? – Giuppox Dec 09 '20 at 07:52
  • The path generation part looks fine. I am not convinced by `calcMSD`. What exactly are you trying to calculate? You are calculating it from a single path, don't you want to average over a bunch of paths? I just do not know the objective here, would be best if you explained that or posted a link to whatever calculation you are trying to replicate. Also, it is likely you are missing a square root as in `msd[i - 1] = np.sqrt(squareddisplacements.mean())` and possibly a division by 'i'? again really depends on what you are trying to calc – piterbarg Dec 09 '20 at 08:15
  • @piterbarg Thanks for you response. I'm trying to calculate the average squared distance from the starting position at every time interval. So if I had 60 timepoints, I'd first calculate if for time intervals of size 1, then of size 2, then all the way to size 60. At 60 there is only 1 step of that size, so it would just calculate the distance from ending point to the starting point. For example, if p(t) gives you position as a function of time, then what I'm trying to calculate at t = 1 would be: [p(1) - p(0))**2 + (p(2) - p(1))**2 + … + (p(n) - p(n - 1))**2] / n – Luke Dec 10 '20 at 04:04
  • @Luke It looks to me that your code does exactly that. Your code can be simplified a bit. Eg your path loop can be simply `for i in range(1, path.size): path[i] = path[i-1] + 2*np.random.randint(0, 2)-1`. Also in `calcMSD` function, `data[:data.size - i]` can be written in a more Pythonic way as `data[:- i]`. otherwise looks good to me – piterbarg Dec 10 '20 at 08:34

1 Answers1

1

Mean Squared Displacement

Definition of mean squared displacement

To make sure we are agreeing on the definition, the mean squared displacement (MSD) is the mean of the squared distance from origin of a collection of particles (in our case, walkers) at a specific step.

Since we are in 1D, the distance is simply the absolute position : distance at time t = | path[t] |

To get the MSD at step t, we must take all the squared positions at step t, then take the mean of the resulting array.

That means : MSD at step t == (path[:, t]**2).mean() (the element-wise squaring makes it useless to take the absolute values)

path[:, t] means here : the positions of each walker at step t.

Generating paths

We must edit your code a bit to have several independant paths.

walkers = 50 # Example for 50 walkers
path = []
steps = 60
for walker in range(walkers):
    walker_path=[0]
    position = 0
    for i in range(1, steps):
        move = np.random.randint(0, 2)
        if move == 0: position += -1
        else: position += 1
        walker_path.append(position)
    path.append(walker_path)
path = np.array(path)

As a sidenote, you could also use random.choice to randomly pick -1 or 1 directly, like so:

for i in range(1, steps):
    position += np.random.choice([-1, 1])
    walker_path.append(position)

From now, path is a 2D array : a row per walker path.

Getting MSD at each step

To get all the MSD at every steps, we vary t from 1 to the number of steps.

You can do it inside a list comprehension in Python (and getting rid of your calcMSD function):

msd = [(path[:, i]**2).mean() for i in range(1, steps)] 
# you could also put range(1, len(path[0]))

Result

Your function was not adapted to plot several walkers.

The proposed solution will give you something like this :

MSD Proposed

What you want

I hope I understood what you truly wanted. I don't know how it's called but it can also be calculated in a list comprehension. Considering the path for a single walker :

not_msd = [(np.array([path[start:start+sep+1] for start in range(steps-sep)])**2).mean() for sep in range(1, steps)]

We go from start to start + step to get the difference between two positions that are step steps apart. We reproduce this operation, starting from 0 to the end of the list. Then we square the obtained list and take the mean.

We do that for all step size : 1 to the whole length.

Result :

Not MSD proposed

Whole Brain
  • 2,097
  • 2
  • 8
  • 18
  • Thanks for a really clear response. I think I've misunderstood how to calculate the MSD. Say p(t) was position as a function of time. If I am understanding you correctly, at t = 3 you would calculate: [ p(1)**2 + p(2)**2 + p(3)**2 ] / 3 I'll try to explain what I was attempting. I thought that you needed to calculate the difference in position between timepoints separated by some number t. So for t = 3 I thought it would be: [ (p(3) - p(0))**2 + (p(4) - p(1))**2 + (p(5) - p(2))**2 + … (p(n) - p(n - 3))**2 ] / (the number of timepoints with 3 steps in between them. – Luke Dec 10 '20 at 04:04
  • My definition of MSD was **wrong**. I confused the mean of each step with the mean of all positions at a specific step. You need several walkers to have a MSD, else the MSD will simply be the square of your single path. I edited my answer. – Whole Brain Dec 10 '20 at 09:13