poisson_point_process.py

import scipy

def PoissonPP( rt, Dx, Dy=None ):
  '''
  Determines the number of events `N` for a rectangular region,
  given the rate `rt` and the dimensions, `Dx`, `Dy`.
  Returns a <2xN> NumPy array.
  '''
  if Dy == None:
    Dy = Dx
    
  N = scipy.stats.poisson( rt*Dx*Dy ).rvs()
  x = scipy.stats.uniform.rvs(0,Dx,((N,1)))
  y = scipy.stats.uniform.rvs(0,Dy,((N,1)))
  P = np.hstack((x,y))
  
  return P


def main():
  rate, Dx = 0.2, 20
  P = PoissonPP( rate, Dx ).T
  fig, ax = subplots()
  ax = fig.add_subplot(111)
  ax.scatter( P[0], P[1], edgecolor='b', facecolor='none', alpha=0.5 )
  # lengths of the axes are functions of `Dx`
  xlim(0,Dx) ; ylim(0,Dx)
  # label the axes and force a 1:1 aspect ratio
  xlabel('X') ; ylabel('Y') ; ax.set_aspect(1)
  title('Poisson Process')
  savefig( 'poisson_lambda_0p2.png', fmt='png', dpi=100 )


if __name__ == '__main__':
  main()