Poisson Process in R from exponential distribution

Question

I am trying to simulate a poisson process sample path in R by starting off with exponentially distributed random variables. For example, for a value of $\lambda=0.5$, I can generate 500 samples and then I want to plot the poisson process path on a time interval of [0,10] for example, how can I do this in R?

attempt

set.seed(1)
n <- 100
x <- cumsum(rexp(50, rate=0.5))
y <- cumsum(c(0, rep(1, 50)))
plot(stepfun(x, y), xlim = c(0, 10), do.points = FALSE, 
     main="L=0.5")

This seems to work, although I don't know how efficient it is since I'm only getting what I want by restricting the x axis of the graph

Answer 1

Basically, you need to compute the successive arrivals $S_i$ for $i=1$, $2$, $\dots$ as cumulative sums of independent exponential interarrivals. So the two main ingredients here are rexp and cumsum. Then you plot the points $[S_i,\, i]$ with a step interpolation (type = "s" in plot functions), and an extra point for $i=0$ and $S_i:=0$ will help.

On a given interval you don't know by advance how many arrivals $S_i$ will come. So you can either use a loop with a break control statement, or simulate more than needed as shown here. The second option may be more efficient in R.

lambda <- 0.5
tMax <- 100

## find the number 'n' of exponential r.vs required by imposing that
## Pr{N(t) <= n} <= 1 - eps for a small 'eps'
n <- qpois(1 - 1e-8, lambda = lambda * tMax)

## simulate exponential interarrivals the
X <- rexp(n = n, rate = lambda)
S <- c(0, cumsum(X))
plot(x = S, y = 0:n, type = "s", xlim = c(0, tMax)) 

## several paths?
nSamp <- 50
## simulate exponential interarrivals
X <- matrix(rexp(n * nSamp, rate = lambda), ncol = nSamp,
            dimnames = list(paste("S", 1:n, sep = ""), paste("samp", 1:nSamp)))
## compute arrivals, and add a fictive arrival 'T0' for t = 0
S <- apply(X, 2, cumsum)
S <- rbind("T0" = rep(0, nSamp), S)
head(S)
## plot using steps
matplot(x = S, y = 0:n, type = "s", col = "darkgray",
        xlim = c(0, tMax),
        main = "Homogeneous Poisson Process paths", xlab = "t", ylab = "N(t)")