def walk():
     'starting at x=0, step left/right with probability 1/2 until x=1'
     x=0
     while x!=1:
         x+=random.choice([-1,1]) # equi-probable left-right move
         yield x

from __future__ import  division # want floating point division

fig,ax=subplots()
p = linspace(0,0.5,20)
ax.plot(p,p/(1-p),label=r'$p<1/2$',lw=3.)
ax.plot([0.5,1],[1,1],label=r'$p > 1/2$',lw=3.)
ax.axis(ymax=1.1)
ax.legend(loc=0)
ax.set_xlabel('$p$',fontsize=16)
ax.set_ylabel('$P$',fontsize=16)
ax.set_title(r'Probability of reaching $x=1$ from $x=0$',fontsize=16)
ax.grid()

def walk():
     'starting at x=0, step left/right with probability 1/2 until x=1'
     x=0
     while x!=1:
         x+=random.choice([-1,1]) # equi-probable left-right move
         yield x

random.seed(123) #  set seed for reproducibility

s = list(walk()) # generate the random steps

fig,ax=subplots()
ax.plot(s)
ax.set_ylabel("particle's x-position",fontsize=16)
ax.set_xlabel('step index k',fontsize=16)
ax.set_title('Example of random walk',fontsize=16)

s = [list(walk()) for i in range(50)] # generate 50 random walks
len_walk=map(len,s) # lengths of each walk

fig,ax=subplots()
for i in s:
    ax.plot(i)
ax.set_ylabel("particle's x-position",fontsize=16)
ax.set_xlabel('step index k',fontsize=16)
ax.set_title('average length=%3.2f'%(mean(len_walk)))

def walk(limit=50):
     'limited version of random walker'
     x=0
     while x!=1 and abs(x)<limit:
         x+=random.choice([-1,1])
         yield x

def nwalk(limit=500):
     'limited version of random walker. Only returns length of path, not path itself'
     n=x=0
     while x!=1 and n<limit:
         x+=random.choice([-1,1])
         n+=1
     return n # return length of walk

len_walk = [nwalk() for i in range(500)] # generate 500 limited random walks

import pandas as pd
from collections import Counter, OrderedDict

lw = pd.Series(Counter(len_walk))/len(len_walk)
fig,ax=subplots()

ax.plot(lw.index,lw.cumsum(),'-o',alpha=0.3)
ax.axis(xmin=-10,ymin=0)
ax.set_title('Estimated Cumulative Distribution',fontsize=16)
ax.grid()

def estimate_std(limit=10,ncount=50):
    'quick estimate of the standard deviation of the averages'
    ws= array([[ nwalk(limit) for i in range(ncount)] for k in range(ncount)])
    return (limit,ws.mean(), ws.mean(1).std())

for limit in [10,20,50,100,200,300,500,1000]:
    print 'limit=%d,\t average = %3.2f,\t std=%3.2f'% estimate_std(limit)

import networkx as nx
import itertools as it

class Graph(nx.DiGraph):
   '''
   operations assuming `pos` attribute in nodes to support drawing and
   manipulating path lattice.
   '''
   def draw(self, ax=None,**kwds):
      '''
      Draw based on `pos` attribute and pass kwds to nx.draw
      :param: axes(optional, default is None)
      '''
      pos = self.get_pos()
      node_size=kwds.pop('node_size',200)
      alpha=kwds.pop('alpha',0.3)
      if ax is None: nx.draw(self,pos=pos,
                             node_size=node_size,
                             alpha=alpha,
                             **kwds)
      else:          nx.draw(self,pos=pos,
                             node_size=node_size,
                             alpha=alpha,
                             ax=ax,**kwds)

   def get_pos(self,n=None):
      '''
      n := str name of node
      get positions as returned dictionary
      '''
      pos=dict([( i,j['pos'] ) for i,j in self.nodes(data=True)])
      if n is None:
         return pos
      else:
         return pos[n]

   def getx(self,x):
      return sorted([(i[0],i[1]['val'] ) for i in self.nodes(True) if i[0][0]==x])

   def gety(self,y):
      return sorted([(i[0],i[1]['val'] ) for i in self.nodes(True) if i[0][1]==y])

# functions to allow diagonal lattice walking
def diagwalk(level,n):
   x = level
   y = -level
   while y<=1 and x<n+1:
      yield (x,y)
      x+=1
      y+=1

def diagwalker(n):
   'daisy-chains the individual diagonal walkers'
   assert n%2 # odd only
   return it.chain(*(diagwalk(i,n) for i in range(n+1)))


def construct_graph(level=5):
    g=Graph()
    g.level = level
    g.add_nodes_from([(i,dict(pos=i)) for i in diagwalker(g.level)])
    
    for x,y in g.nodes():
       if y!=1 and x< g.level:
          g.add_edge((x,y), (x+1,y-1) )
          g.add_edge((x,y), (x+1,y+1) )
    
    g.node[(0,0)]['val']=1L # long int
    for j in diagwalker(g.level):
        if j==(0,0): continue
        x,y=j
        g.node[j]['val']=sum(g.node[k]['val'] for k in g.predecessors(j))
    return g
g=construct_graph()

fig,ax=subplots()
fig.set_size_inches((6,6))
ax.set_aspect(1)
ax.set_title('Directed Path Lattice',fontsize=18)
g.draw(ax=ax)

l5=list(nx.all_simple_paths(g,(0,0),(5,1)))
for i in l5:
    print i

fig,ax=subplots()
fig.set_size_inches((6,6))
g.draw(ax,with_labels=0)
g.subgraph(l5[0]).draw(ax=ax,with_labels=1,node_color='b',node_size=700)
g.subgraph(l5[1]).draw(ax=ax,with_labels=1,node_color='b',node_size=800)
ax.set_title('Pathways that terminate in five steps',fontsize=16)

for i in g.getx(5):
    print 'node (%d,%d)\t'%(i[0]),
    print 'number of paths = %d'%(i[1])

import operator as op

def get_cond_prob(g):
   'compute conditional probability for later' 
   cp = OrderedDict()
   for i,j in g.gety(1):
      if (i[0],-i[0] ) not in g: continue
      cp[i[0]] = g.node[i]['val']/sum(map(op.itemgetter(1),g.getx(i[0])))
   return cp

def get_prob(g,cp=None):
   'get unconditional probability of termination '
   prob = OrderedDict()
   cq = OrderedDict()
   if cp is None: cp = get_cond_prob(g)
   for i,j in cp.iteritems(): cq[i]=1-j
   prob[1]=  0.5  # initial condition
   for i in cp.keys():
      prob[i]= cp[i]*prod(cq.values()[:(i-1)//2][::-1])
   return prob


g=construct_graph(15)
p=get_prob(g)
print 'prob of termination in 5 steps = ',p[5]

fig,ax=subplots()
ax.axis(ymax=1)
ax.plot(p.keys(),p.values(),'-o')
ax.set_xlabel('exact no. of steps to terminate',fontsize=18)
ax.set_ylabel('probability',fontsize=18)
ax.axis(ymax=.6)

import sympy

# analytical result from Feller
pfeller = dict([(2*k-1,sympy.binomial(1/2,k)*(-1)**(k-1)) for k in range(1,10)])

for k,v in p.iteritems():
    assert v==pfeller[k] # matches every item!

def nwalk(limit=500):
     'p=2/3 version of random walker. Only returns length of path'
     n=x=0
     while x!=1 and n<limit:
         x+=random.choice([-1,1,1]) # twice as many 1's as before
         n+=1
     return n # return length of walk

mean([nwalk() for i in range(100)])

for limit in [10,20,50,100,200,300,500,1000]:
    print 'limit=%d,\t average = %3.2f,\t std=%3.2f'% estimate_std(limit,100)