#test the RWG simulator
nTrials = 100000
data = [randint(0,1) for _ in range(nTrials)]  #a list of 100,000 random zeros (tails) and ones (heads)
print(data.count(0), data.count(1)) # how many of each? Should be about equal...

#statistics of flipping a fair coin
def statsCoin(nFlips):
    flips = [randint(0,1) for _ in range(nFlips)]
    first = flips[0]
    last = flips[-1]
    heads = flips.count(1)
    #compute longest streak
    streak = 1 #first flip starts a streak
    biggestStreak = streak # so far so good
    start = flips[0]
    pos = 1
    while pos < nFlips:
        if flips[pos]==start:
            streak += 1
            if streak > biggestStreak:
                biggestStreak = streak
        else:
            streak = 1
        start = flips[pos]
        pos += 1
    return (first, last,heads, biggestStreak)

nTrials = 100000
data=[]
for k in range(nTrials):
    data.append((statsCoin(10)))
firstData=[d[0] for d in data]
headsData=[d[2] for d in data]
streakData=[d[3] for d in data]
[float(streakData.count(i)/nTrials) for i in [1..10]]

[float(firstData.count(i)/nTrials) for i in range(1)]

[float(headsData.count(i)/nTrials) for i in range(10)]

sum([0.00096,
 0.00955,
 0.04364,
 0.11691,
 0.20541,
 0.24755,
 0.20348,
 0.1166,
 0.04449,
 0.01031])

3

#fair random walk simulation
# how long to first return?
# assume we start t=0 at location =0

def time_to_first_return():
    location = 2*randint(0,1)-1
    t = 1
    while  not location==0:
        t += 1
        location += 2*randint(0,1)-1
    return t

nTrials = 1000
data = [time_to_first_return() for _ in range(nTrials)]
print(data)

[data.count(2*k) for k in [1..20]]

#fair random walk simulation
# how long to first return?
# assume we start t=0 at location =0

def time_to_first_return():
    location = 2*randint(0,1)-1
    t = 1
    while  not location==0:
        t += 1
        location += 2*randint(0,1)-1
    return t

nTrials = 200
data = [time_to_first_return() for _ in range(nTrials)]
print(data)

[data.count(2*k) for k in [1..20]]

big=[a for a in data if a>1000]
big.sort()
big

# simple 1-dim RW
nSteps = 10000
nWalks = 8
canvas = Graphics()

for w in range(nWalks):
    y = 0
    path=[y]
    for _ in range(nSteps):
        y += 2*randint(0,1)-1
        path.append(y)
    canvas += list_plot(path,plotjoined=True, color=Color(w/nWalks,1-w/nWalks,.2))
show(canvas)

#gather lotsa data
nSteps = 1000
nTrials = 1000
#collect end distance, number of returns, first return
end_dist =[]
num_ret = []
first_ret = [] 

for _ in range(nTrials):
    end=0; nr=0; fr=0 #fr will change from 0 if there is a return
    for step in [1..nSteps]:
        end += 2*randint(0,1)-1
        if end == 0:
            nr += 1
            if nr ==1:
                fr = step
    end_dist.append(abs(end))
    num_ret.append(nr)
    first_ret.append(fr)
    

for k in [10..20]:
    print (end_dist[k], num_ret[k], first_ret[k])

histogram(end_dist)

hret1000=histogram(num_ret)

sum(end_dist)

hret1000.save("histret1000.pdf")

sum(end_dist)

histogram(num_ret)

h1000ret=histogram(num_ret)
h1000ret.save("histo1000numret.pdf")

num_ret.count(0)

histogram(first_ret)

first_ret.count(0)

max(first_ret)

import numpy as np
nSteps = 100
nTrials = 1000
def b(n):
    return float(binomial(2*n,n)/2^(2*n))
#collect end distance, number of returns, first return
end_dist =[]
#num_ret = []
#first_ret = [] 

for _ in range(nTrials):
    end=0#fr will change from 0 if there is a return
    for step in [1..nSteps]:
        end += 2*randint(0,1)-1
        
    end_dist.append(abs(end))
print(np.mean(end_dist),b(nSteps/2)*nSteps)