1

我正在模拟二维随机游走，方向 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？

4

2 回答 2

4

这是计算均方位移 (MSD) 的快速片段。路径由时间间隔相等的点组成，就像您的 randwalk 的情况一样。您可以只放置 12-walk 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
于 2014-12-30T16:46:59.923 回答
0

首先，您实际上不需要存储整个 1000 步步行，只需要存储最终位置。

此外，没有理由将它们存储到文本文件并重新加载它们，您可以在内存中使用它们 - 只需将它们放在数组列表中，或者放在多一维的数组中。即使您需要将它们写出来，您也可以这样做并保留最终值，而不是代替. （此外，如果您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)
于 2014-10-20T18:46:03.390 回答