python,scipy- 作业 -- 如何在网格中保留粒子位置历史

Question

我遇到了以下问题：

我要创建一个程序，让一个粒子在两种情况下进行随机行走，总共走N=1000步：一是在线性的一维空间中，二是在二维空间中。这个程序需要计算一个平均值(S)，其中S是粒子至少访问过一次的网格位置的数量。我要进行10000次实验，并找出10个点（每100步一个，从0到1000），这些点将是10000次实验的平均值。最后，绘制平均值(S)与时间t的关系图。

我写了以下代码：

import scipy as sc
import matplotlib.pyplot as plt
import random

plegma=1000
grid=sc.ones(plegma)   # grid full of available positions(ones)    

for p in range(10000):
    #-------------------Initialize problem-------------------------------------------------
    his_pos=[]                  # list which holds the position of the particle in the grid
    in_pos = int(sc.random.randint(0,len(grid),1))            #initial position of particle
    means=[]                                                    #list which holds the means 
    #--------------------------------------------------------------------------------------
    for i in range(0,1000,100):
        step=2*sc.random.random_integers(0,1)-1        #the step of the particle can be -1 or 1
        # Check position for edges and fix if required
        # Move by step
        in_pos += step
        #Correct according to periodic boundaries
        in_pos = in_pos % len(grid)  

        #Keep track of random walk
        his_pos.append(in_pos)
        history=sc.array(his_pos)
        mean_his=sc.mean(history) 
        means.append(mean_his)


plt.plot(means,'bo')
plt.show()

更新 -------------------------------------

import scipy as sc
import matplotlib.pyplot as plt
import random

plegma=1000
his_pos=[] # list which holds the number of visited cells in the grid
means=[] #list which holds the means

for p in range(10000):
    #-------------------Initialize problem-------------------------------------------------
    grid=sc.ones(plegma)   # grid full of available positions(ones)      
    in_pos = int(sc.random.randint(0,len(grid),1))            #initial position of particle
    num_cells=[]       # list which holds number of visited cells during run                         
    #--------------------------------------------------------------------------------------
    for i in range(1000):
        step=2*sc.random.random_integers(0,1)-1 #the step of the particle can be -1 or 1

        # Check position for edges and fix if required
        # Move by step
        in_pos += step
        #Correct according to periodic boundaries
        in_pos = in_pos % len(grid)  

        grid[in_pos]=0  # mark particle position on grid as "visited"

        if (i+1) % 100 == 0:

            number=1000-sc.sum(grid)  # count the number of "visited positions" in grid
            num_cells.append(number)  # append it to num_cells

      his_pos.append(num_cells) # append num_cells to his_pos
      history=sc.array(his_pos)


mean_his=history.mean(1)
means.append(mean_his)

更新2 -----------------------------

    if (i+1) % 10 == 0:

            number=1000-sc.sum(grid)  # count the number of "visited positions" in grid
            num_cells.append(number)  # append it to num_cells

    his_pos.append(num_cells) # append num_cells to his_pos
    history=sc.array(his_pos)
mean_his=history.mean(0)

plt.plot(mean_his,'bo')
plt.show()

谢谢！

Answer 1

我不太明白问题具体在问什么（时间 t 在你的公式中是怎么考虑的？point 指的是什么？），不过关于你想进行的操作，我明白你需要把每次迭代得到的10个平均值数组 mean_his 都放到一个最终的10000x10的 means 数组里。

每个 mean_his 数组是从一个包含100个步骤的一维数组中准备的。我用一个包含十个步骤的数组来举个例子，这个数组每两个步骤取一次平均（而不是每1000个步骤取一次平均）：

>>> his_pos = [1,2,3,4,5,6,7,8,9,10]    #the ten positions
>>> history = np.array(his_pos)
>>> history
array([ 1,  2,  3,  4,  5,  6,  7,  8,  9, 10])
>>> ar2 = np.reshape(history, (-1,2))   # group in two in order to calculate mean
>>> ar2
array([[ 1,  2],
       [ 3,  4],
       [ 5,  6],
       [ 7,  8],
       [ 9, 10]])
>>> mean_his = ar2.mean(1)
>>> mean_his
array([ 1.5,  3.5,  5.5,  7.5,  9.5])
>>>  

然后你需要把 mean_his 添加到 means 中10000次，并且以类似的方式计算平均值（注意，means 需要在外层循环外初始化，这样就不会在每次重复时重新初始化）。

Answer 2

1) 是的，10000步是指你想要的准确度——你需要在时间点 t = 100, 200 ... 1000 时，计算10000次随机走动中访问过的格子的平均数量。为了得到这些数据，你需要记录每次随机走动中访问过的格子数量（总共10000次）。为此，你需要在 p 循环外先初始化你的问题（也就是初始化 his_pos 和 means）。在 p 循环内部，你应该初始化你的随机走动——获取你的网格、初始位置和你要写入结果的列表。所以，问题的初始化大概是这样的：

import scipy as sc
import matplotlib.pyplot as plt

plegma=1000
his_pos=[]  # list which holds the position of the particle in the grid
means=[]    #list which holds the means 

for p in range(10000):
    #-------Initialize problem--------
    in_pos = int(sc.random.randint(0,len(grid),1)) #initial position of particle
    grid=sc.ones(plegma)   # grid full of available positions(ones)    
    his_pos.append([])

初始化后，我们需要进行随机走动，也就是 i 循环。目前你只进行了10步随机走动（len(range(0,1000,100)) == 10），但实际上走动应该有1000步。因此，i 循环应该是：

    for i in range(1000):

在走动过程中，你需要以某种方式标记访问过的格子。最简单的方法是把 grid[in_pos] 改成0。然后，每走到第100步时，你需要计算一下访问过的格子数量。实现这个的方式是：

        if i % 100 == 0:
            # count the number of 0s in grid and append it to his_pos[p]

最后，在你完成1000次随机走动后，你的 his_pos 将会是一个包含(10000*10)的列表的列表，你需要从中按列计算平均值。

更新：

为了在多次运行中不丢失信息，我们应该把存储 p 次运行结果的列表，添加到所有结果的列表中。后者就是 his_pos。为了实现这一点，我们可以在 his_pos 中添加一个空列表，并在 p 次运行中填充结果，或者在 p 次运行之前初始化一个空列表，并在 p 次运行之后将其添加到 his_pos 中。算法大概会是这样的：

#-------Initialize problem--------
plegma=1000
his_pos=[]  # list which holds the number of visited cells in the grid
means=[]    #list which holds the means 

for p in range(10000):
    #-------Initialize run--------
    in_pos = int(sc.random.randint(0,len(grid),1)) #initial position of particle
    grid=sc.ones(plegma)   # grid full of available positions(ones)    
    num_cells = []  # list which holds number of visited cells during run  
    for i in range(1000):
        # make i-th step, get particle position
        # mark particle position on grid as "visited"
        if (i+1) % 100 == 0: # on each 100th step (steps count from 0, thus i+1)
            # count the number of "visited positions" in grid
            # append it to num_cells
    # append num_cells to his_pos
# find column-wise means for his_pos