Вычисление среднеквадратичного смещения двумерного случайного блуждания в Python

Question

Я моделирую двухмерное случайное блуждание с направлением 0 <θ <2π и T = 1000 шагов. У меня уже есть код, который имитирует одну прогулку, повторяет ее 12 раз и сохраняет каждый прогон в последовательно именованные текстовые файлы:

a=np.zeros((1000,2), dtype=np.float)
print a                                   # Prints array with zeros as entries

# Single random walk
def randwalk(x,y):              # Defines the randwalk function
    theta=2*math.pi*rd.rand() 
    x+=math.cos(theta);
    y+=math.sin(theta);
    return (x,y)                # Function returns new (x,y) coordinates

x, y = 0., 0.                   # Starting point is the origin
for i in range(1000):           # Walk contains 1000 steps
    x, y = randwalk(x,y)
    a[i,:] = x, y               # Replaces entries of a with (x,y) coordinates

# Repeating random walk 12 times
fn_base = "random_walk_%i.txt"      # Saves each run to sequentially named .txt
for j in range(12):
    rd.seed()                       # Uses different random seed for every run
    x, y = 0., 0.
    for i in range(1000):
        x, y = randwalk(x,y)
        a[i,:] = x, y
    fn = fn_base % j                # Allocates fn to the numbered file
    np.savetxt(fn, a)               # Saves run data to appropriate text file

Теперь я хочу рассчитать среднеквадратичное смещение за все 12 прогулок. Для этого моей первоначальной мыслью было импортировать данные из каждого текстового файла обратно в массив Numpy, например:

infile="random_walk_0.txt"
rw0dat=np.genfromtxt(infile)
print rw0dat

А затем каким-то образом манипулировать массивами, чтобы найти среднеквадратичное смещение.

Есть ли более эффективный способ найти MSD с тем, что у меня есть?

Answer 1

Вот быстрый фрагмент для вычисления среднеквадратичного смещения (MSD). Там, где путь состоит из точек, равномерно распределенных во времени, как, кажется, имеет место для вашего рандома. Вы можете просто поместить в 12-шаговый цикл for и вычислить его для каждого a[i,:]

#input path =[ [x1,y1], ... ,[xn,yn] ].

def compute_MSD(path):
   totalsize=len(path)
   msd=[]
   for i in range(totalsize-1):
       j=i+1
       msd.append(np.sum((path[0:-j]-path[j::])**2)/float(totalsize-j))

   msd=np.array(msd)
   return msd

Answer 2

Во-первых, вам на самом деле не нужно хранить всю 1000-шаговую прогулку, а только конечную позицию.

Кроме того, нет причин хранить их в текстовых файлах и загружать их обратно, вы можете просто использовать их в памяти - просто поместите их в список массивов или в массив еще на 1 измерение. Даже если вам нужно выписать их, вы можете сделать это, а также сохранить окончательные значения вместо. (Кроме того, если вы на самом деле не используете numpy для производительности или простоты построения 2D-массива вы можете рассмотреть возможность его построения итеративно, например, используя csv модуль, но это еще один вызов.)

В любом случае, учитывая ваши 12 последних позиций, вы просто рассчитываете расстояние от каждой из (0, 0) затем возведите в квадрат, суммируйте их все и разделите на 12. (Или, поскольку очевидный способ вычислить расстояние от (0, 0) это просто добавить квадраты x а также y позиции и затем квадратный корень результат, просто пропустите квадратный корень и квадрат в конце.)

Но если по какой-то причине вы хотите сохранить каждую прогулку в файле, то после загрузки их обратно, walk[-1] дает вам окончательную позицию в виде одномерного массива из 2 значений. Таким образом, вы можете либо прочитать эти 12 конечных позиций в массив 12x2 и векторизовать среднеквадратичное расстояние, либо просто накапливать их в списке и делать это вручную.

Пока мы на это, rd.seed() не обязательно; весь смысл PRNG заключается в том, что вы продолжаете получать разные числа, если вы явно не сбросите начальное значение к его первоначальному значению, чтобы повторить их.

Вот пример отбрасывания двух дополнительных сложностей и выполнения всего напрямую:

destinations = np.zeros((12, 2), dtype=np.float)
for j in range(12):
    x, y = 0., 0.
    for i in range(1000):
        x, y = randwalk(x, y)
    destinations[j] = x, y

square_distances = destinations[:,0] ** 2 + destinations[:,1] ** 2
mean_square_distance = np.mean(square_distances)