from scipy import *
from pylab import *
from scipy.misc import *

NumSteps = 20
Seed = 10000
seed(Seed)
NumAve = 10000
Position = zeros((NumAve))

cla()

bins = array([i for i in range(-NumSteps,NumSteps)])

#1. Create NumAve random walks and record their final endpoint
for i in range(NumAve):
    N_rand = randint(0,2, NumSteps)
    dir_rand = 2*N_rand -1
    Position[i] = sum(dir_rand)

#2. Make a histogram of the endpoints
pdf, bins, patches = hist(Position, bins, normed=1, rwidth=1, align='left')

#if you don't understand this line, type it into ipython, and then type "y" on a new line
y = array([i for i in range(-NumSteps/2,NumSteps/2)])

#3. Compare the result with the exact result obtained by applying the binomial distribution
prefactor   = factorial(NumSteps)*(0.5)**NumSteps
distr # = prefactor/?

#4. Plot the results
plot(2*y,distr)

#5. Compare the exact result to that obtained for a gaussian with the correct variance
gauss # = depends on NumSteps and y

plot(2*y,gauss)


show()