# poisson-process.R
# Bret Larget
# January 23, 2009

# For Statistics 310

# One characterization of the Poisson process
# is that inter-event times are iid Exponential random variables.
# Here is code to simulate a Poisson process and to count the
# number of events before a given time.  These counts
# can then be compared to the exact Poisson distribution.

################################################################################
# getN() --- function to count the number of events before a specified time
# arguments:
#   x = a vector of inter-event times (assumed to be positive)
#   tt = a positive time
# return value:
#   the number of events before time tt (possibly 0)
################################################################################

getN = function(x,tt=1) {
  if(sum(x)<tt)
    return(NA)
  if(x[1]>tt)
    return(0)
  else
    return(max(which(cumsum(x)<tt)))
}

# How the code works:
#  If the sum of all of the x values is less than tt,
#  getN() returns NA for a missing value
#  x[1] is the first element of x.
#  In general, the [ ] operator can be used in powerful ways to select subsets of an object
#  The line
#         return(max(which(cumsum(x)<tt)))
#  deserves some explanation.  Start at the inside and work out.

#  1. The expression cumsum(x) is a vector the same length of x where
#     each element is the cumulative sum.
#  2. (cumsum(x) < tt) is a vector where each element is replaced by T
#     or F depending on whether or not the expression is true or false.
#  3. which() returns the indices of a vector that evaluate as true; the input must be logical
#     Here, if the input x is all positive, all of the T elements will be at the beginning.
#  4. max() returns the maximum; here, the maximum corresponds to the last event before tt

# EXAMPLE:
# Create an initial vector of time intervals
# Notice that the event times are 0.1, 0.7, 1.2, 1.5, 2.9
# > x = c(0.1,0.6,0.5,0.3,1.4)
#
# Call getN() on this vector; by default, the time is 1
# In the output, the [1] indicates that the row begins with the first element
# There are exactly two event times before the time 1
# > getN(x)
# [1] 2
#
# Now call getN() with the second argument equal to 2
# Notice that there are four events before time 2
# > getN(x,2)
# [1] 4
#
# Now call getN() for time 5.
# Since all events occur before time 5, we do not know of the next event would be before or after time 5.
# Hence, the program returns NA for a missing value.
# > getN(x,5)
# [1] NA
#
# The cumsum() function computes the cumulative sum of a vector; here the event times.
# > cumsum(x)
# [1] 0.1 0.7 1.2 1.5 2.9

################################################################################
# simPP() --- function to simulate a Poisson process many times and return counts
# arguments:
#   ntrials = the number of times to repeat the Poisson process simulation
#   tt = the length of the time interval
#   lambda = the rate parameter
# return value:
#   a vector of counts, one for each trial, of the number of events before time tt
#
################################################################################

simPP = function(ntrials,tt=1,lambda=2) {
  maxN = round(qpois(0.99999,tt*lambda))+10 # large enough so that the chance of needing a larger sample is tiny
  x = matrix(rexp(ntrials*maxN,lambda),ntrials,maxN)
  return(apply(x,1,getN,tt=tt))
}

# How the code works:
# The general idea is to create a matrix of iid random exponential random variables
# so that each row corresponds to a single realization of a Poisson process.
# Enough random variable are generated so that it is very improbable that the sum in any row will be below the target time tt.
# So, x is an ntrials by maxN matrix.
# The function apply() is used to apply a function to either the rows or columns of a matrix.
# The command apply(x,1,getN,tt=tt) instructs R to take the matrix x
# and apply to each row (1 stands for the first dimension of the two-dimensional array, (row,column))
# the function getN.  The additional argument tt=tt is passed on to getN().
# Here, the left tt matches the name of an argument in getN() and the right tt matches the name
# of an object already defined in simPP(), here as one of its arguments.
# So, this code applies getN() to each row of x and returns the answer as a vector of length ntrials,
# the number of rows of x.