Simulation Poisson Process using R and ggplot2

Question

Using a simulation of a Poisson process with rate lambda = 0.7. Show a sample run of a Poisson process with N(t) on the vertical axis and time t on the horizontal axis. The simulation is from the range t[0:100]. Generate first graph with 10 trajectories and second graph with 100 trajectories.

I have tried the following code but I cannot generate both graphs.

library(plyr)
library(ggplot2)

Process_poisson<- function(t, lambda){
distr_poisson<- rpois(1, t*lambda)
s_poisson<- sort(runif(distr_poisson, 0, t))
data.frame(x = c(0, 0, s_poisson),y = c(0, 0:distr_poisson))
}

N_simulations<- function(n,t,lambda){
s_poisson<- lapply (1:n, function(n) data.frame(Process_poisson(t, lambda), simulation = n))
s_poisson<- ldply (s_poisson, data.frame)
s_poisson$simulation<- factor(s_poisson$simulation)
}

t<- 0:100
lambda<- 0.7
N_simulations(10, t, lambda)
N_simulations(100, t, lambda)

par(mfrow = c(1,2))

matplot(x, y, type = "l", lty = 0:5, lwd = 1, lend = par("lend"),
     pch = NULL, col = simulation, cex = 0.5, bg = NA, main =sprintf("Nº simulations of trajectories of Poisson Process",10,lambda), xlab = "Time", ylab = "N(t)",
   xlim = c(0,100), ylim = c(-10,0))

matplot(Proceso_poisson(t, lambda), n, y, type = "l", lty = 0:5, lwd = 1, lend = par("lend"),
     pch = NULL, col = simulacion, cex = 0.5, bg = NA, main =sprintf("Nº simulations of trajectories of Poisson Process",10,lambda), xlab = "Time", ylab = "N(t)",
     xlim = c(0,100), ylim = c(-10,0))

How could I do it?

Thanks so much!

Answer 1

I think you could make this simpler. Here's a ggplot solution.

First, create a function that will simulate a Poisson process by taking samples drawn from an exponential distribution with the appropriate lambda. In this example I have used a while loop that starts with a vector x whose first element is 0. The function grows this vector by adding random samples until its sum reaches the target duration tmax. This is not the most efficient way to do it, but should make the example clearer.

When the target is reached, the function returns the cumulative sum of the vector, which represents the arrival times of a Poisson process of the appropriate lambda. Note that to make plotting easier, it actually returns a data frame with the cumulative times, the cumulative count, and a grouping variable run that will allow us to plot several runs easily on a single plot.

make_sample_df <- function(run, tmax, lambda)
{
  x <- 0
  while(sum(x) < tmax) x <- c(x, rexp(1, lambda))
  data.frame(t = cumsum(x), N = seq_along(x), run = rep(run, length(x)))
}

We can now use this function inside our actual plotting function:

plot_poisson <- function(runs, tmax, lambda)
{
  # Creates one data frame for each run, this sticks them all together:
  df <- do.call("rbind", lapply(seq(runs), make_sample_df, tmax, lambda))

  ggplot2::ggplot(df, aes(t, N, group = run)) + 
    geom_step(alpha = 0.25) + 
    labs( title = paste(runs, "runs of Poisson process with lambda", lambda)) +
    theme(legend.position = "none") +
    coord_cartesian(xlim = c(0, tmax))
}

So you can do:

plot_poisson(runs = 10, tmax = 100, lambda = 0.7)

plot_poisson(runs = 100, tmax = 100, lambda = 0.7)