4
$\begingroup$

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 = F,main="L=0.5")

This seems to work, although i dont know how efficient it is since im only getting what i want by restricting the x axis of the graph

$\endgroup$
5
  • $\begingroup$ You mean something along the lines of hist(rpois(500, lambda = 0.5))? $\endgroup$ Apr 30 '15 at 6:05
  • 1
    $\begingroup$ @Roman I'm pretty sure that's not what's being asked there. $\endgroup$
    – Glen_b
    Apr 30 '15 at 11:08
  • 1
    $\begingroup$ dimebucker -- when you say "plot the Poisson process path" do you mean some kind of step function that jumps every time there's an event in the process? $\endgroup$
    – Glen_b
    Apr 30 '15 at 11:09
  • $\begingroup$ @Glen_b yes that's what I meant $\endgroup$ Apr 30 '15 at 11:51
  • $\begingroup$ @dimebucker91 I think your x and y's are backwards here. $\endgroup$ Mar 4 '19 at 1:10
6
$\begingroup$

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)")
$\endgroup$
1
  • $\begingroup$ @dimebucker91. Sorry, I did not see that you had given nearly the same answer yourself in the question. $\endgroup$
    – Yves
    Apr 30 '15 at 12:21

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.