2

i have the following problem:

Create a program where a particle will execute a random walk for N=1000 steps for these two cases: i) in a 1D system ii) in a 2D system. The program must calculate the mean(S) where S is the number of grid positions where the particle visited at least one time.You will make 10000 runs and find 10 points (one for every 100 steps ,from 0 to 1000) , which will be the means of 10000 runs.Do the plot of mean(S) in relation to time t.

I did this code:

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()

UPDATED -------------------------------------

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)

UPDATE 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()

Thanks!

George
  • 5,808
  • 15
  • 83
  • 160
  • If this is a homework problem (and it sure looks like one) then please flag it as such. – talonmies Nov 14 '11 at 18:00
  • "The program must calculate the where S is the..." Could you please reedit your text to make it understandable? – joaquin Nov 16 '11 at 20:22
  • @joaquin:Ok,i edit that and i updated the post.If you could explain me sth from my questions i'll be grateful! – George Nov 16 '11 at 20:59
  • Should link http://stackoverflow.com/questions/8038420/python-construction-of-lattice-which-traps-molecules-doesnt-work-right/8040549#8040549 and http://stackoverflow.com/questions/8067800/scipy-construction-of-lattice-which-traps-molecules-in-2d-dimension – Benjamin Nov 16 '11 at 23:03
  • @Benjamin:Hello,this is a different exercise.The above you are writting are linked between them.(And as you can see ,i used your tips :) ) – George Nov 17 '11 at 09:28

2 Answers2

0

I am not sure of understanding what the problem is asking for (where time t is considered in your equations?, what do point refer to?), but relative to the operations you are trying to carry out, I understand you have to include each or your 10-means mean_his arrays resulting from each iteration in a final 10000x10 means array.

Each mean_his array is prepared from an 1D array with the 100 steps. I will exemplify this with a array of ten steps that have to be averaged every two (instead of 1000 every 100):

>>> 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])
>>>  

then you append mean_his to means 10000 times and calculate similarly the mean (Note that means have to be initialized outside of the outter loop in order not to be initialized every repetition).

joaquin
  • 82,968
  • 29
  • 138
  • 152
0

1) Yes, 10000 steps refer to desired accuracy - you have to get the mean number of visited cells at time t = 100, 200 ... 1000 for 10000 random walks. To get the data, you have to accumulate the number of visited cells for each random walk you do (of 10000). to do this, you have to initialize your problem (that is, initialize his_pos and means) outside of p loop. Inside of p loop you should initialize your random walk - get your grid, initial positions and the list you'll write results to. So, problem init will look something like

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([])

After init we need to perform random walk - the i loop. Presently you make just 10 steps of random walk (len(range(0,1000,100)) == 10), but the walk should contain 1000 steps. Hence, the i loop should be

    for i in range(1000):

During the walk you have to mark visited cells somehow. The simplest way is to change grid[in_pos] to 0. Then, every 100th step you need to count the number of visited cells. The way to achieve this is like

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

Finally, at the end of your 1000 random walks your his_pos will be the (10000*10) list of lists, out of which column-wise means should be taken.

Update:

In order for not to lose information throughout runs, we should append the list where results for p-th run are stored, to the list with all results. The latter is his_pos. To achieve this, we can either append empty list to his_pos and populate it with results during p-th run, or initialize empty list before p-th run and append it to his_pos after p-th run. The algorithm then will look as follows:

#-------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
Andrey Sobolev
  • 12,353
  • 3
  • 48
  • 52
  • :Hello and thanks for the tips.I didn't understand the 'his_pos.append([])' in the p loop and how to append the count of zeros at the end.I updated my post,i edited the code according to your advice. – George Nov 18 '11 at 20:43
  • We need to get the average of a value over 10000 runs, right? So, the values must be accumulated somehow. In this example I take `his_pos`, which is the accumulator of results, and append it with an empty list, which will accumulate results for the `p`-th run. Equally, you can initialize this accumulator *before* `p`-th run, and append it to `his_pos` *after* `p`-th run. I'll edit the answer to show this opportunity. As for the count of zeros, first you should mark some cells as zeros (you do not do that even in your revised code). Then you can count zeros simply as `1000 - sc.sum(grid)`. – Andrey Sobolev Nov 19 '11 at 04:37
  • :I updated my code again,i think i am close now..You said about finding column-wise means for his_pos ,but his_pos is 1D. – George Nov 19 '11 at 11:01
  • Yep, almost there. Now you should pay attention to 2 facts: 1) you need to find average values over 10000 runs. That means, that averaging should be done when 10000 runs are over, or *after* `p` loop. 2) you append `his_pos` with list every `p`-th run, so `his_pos` is list of lists, hence 2D. You can check that by printing `his_pos` after `p` loop. When you just calculate `sc.mean(his_pos)`, you get average over whole array. To understand how to get column-wise means, you should check joaquin's answer below. – Andrey Sobolev Nov 19 '11 at 15:58
  • :I did that but it isn't correct.It show a 'line' of points.I want to show me 10 points which will be like a curve.Sth am i doing wrong here? – George Nov 19 '11 at 16:18
  • Ok, did you get 10000*10 `his_pos` list? What does `print his_pos` show you? – Andrey Sobolev Nov 19 '11 at 16:54
  • The his_pos is sth like "[[12.0], [15.0], [22.0], [16.0],...." and it's len is 10000.The shape of 'means' is (1,10000).The shape of 'history 'is (10000,1) – George Nov 19 '11 at 17:13
  • That's because you took row-wise averages. To get column-wise averages, you should change axis in `sc.mean` – Andrey Sobolev Nov 19 '11 at 17:29
  • :I am confused now!I tried 'history=sc.transpose(history)' but still the same! – George Nov 19 '11 at 17:55
  • Strange, because your code (with insertion of the `transpose` line in the right place in code) works for me. Then pls check that `len(num_cells)` is equal to 10 before appending it to `his_pos`. If not, then smth is wrong with the length of random walk (the `i` loop) or with `if` line. – Andrey Sobolev Nov 19 '11 at 18:40
  • :I updated my post!In order to have len(num_cells)=10 ,you must use %10.Now, (and without the transpose) the plot is right (the points draw a curve ) but the x-axis is from 0 to 9.I want the axis to be from 0 to 999 (with the plot only 10 points).I must somehow keep in my calculations an array which will have the number of the steps(1000) inside it! – George Nov 19 '11 at 19:13
  • To get the right plot, you need to supply `x` array to `plt.plot`: `plt.plot(sc.arange(100,1000,100),mean_his,'bo')` – Andrey Sobolev Nov 19 '11 at 19:35
  • :The problem is that i have to correlate this plot with the plot that comes from expression 'sn=(8*steps/sc.pi)**(1/2.0)*(1+1/4*steps-3/64*(steps**2))' and my points must lie in that curve.But my points are a lot away from the curve,so sth we are doing wrong? – George Nov 20 '11 at 16:09
  • Check once again - the length of random walk (the `i` loop) should be 1000 (you probably have 100), and results are to be taken (the `if` line) every 100th step (you for sure have 10). Your updated code works for me, and the `mean_his` points lie on the line with the error of <0.05. – Andrey Sobolev Nov 20 '11 at 16:48
  • :Yes!You we right.I had it to 100!(i don't remember why!).Thank you very much for your help!(if sometime you could check a question here 'http://stackoverflow.com/questions/8098065/python-matplotlib-calculate-survival-probability ' i'll be grateful!!It works fine but i am taking very small values in y-axis!Thanks again! – George Nov 20 '11 at 17:06