Solved – Manually simulating Poisson Process in R

poisson processrstochastic-processes

The following problem tells us to generate a Poisson process step by step from $\rho$ (inter-arrival time), and $\tau$ (arrival time).

One of the theoretical results presented in the lectures gives the
following direct method for simulating Poisson process:

• Let $τ_0 = 0$.
• Generate i.i.d. exponential random variables $ρ_1, ρ_2, . . ..$
• Let $τ_n = ρ_1 + . . . + ρ_n$ for $n = 1, 2, . . . .$
• For each $k = 0, 1, . . .,$ let
$N_t = k$ for $τ_k ≤ t < τ_{k+1}$.

  1. Using this method, generate a realization of a Poisson process $(N_t)_t$ with $λ = 0.5$ on the interval $[0, 20]$.
  2. Generate $10000$ realizations of a Poisson process $(N_t)_t$ with $λ = 0.5$ and use your results to estimate $E(N_t)$ and $Var(N_t)$. Compare the estimates
    with the theoretical values.

My attempted solution:

First, I have generated the values of $\rho$ using rexp() function in R.

rhos <-function(lambda, max1)
{
    vec <- vector()

    for (i in 1:max1) 
    {
        vec[i] <- rexp(0.5)
    }

    return (vec)
}

then, I created $\tau$s by progressive summing of $\rho$s. 

taos <- function(lambda, max)
{
    rho_vec <- rhos(lambda, max)
    #print(rho_vec)

    vec <- vector()
    vec[1] <- 0
    sum <- 0
    for(i in 2:max)
    {
        sum <- sum + rho_vec[i]
        vec[i] <- sum
    }

    return (vec)
}

The following function is for finding the value of $N_t = k$ when the value of $k$ is given. Say, it is $7$, etc.

Ntk <- function(lambda, max, k)
{
    tao_vec <- taos(lambda, max)
    val <- max(tao_vec[tao_vec < k])
}

y <- taos(0.5, 20)
x <- seq(0, 20-1, by=1)

plot(x,y, type="s")

Output:

enter image description here

As you can see, the plot of the Poisson process is blank rather than a staircase.

If I change rexp to exp, I get the following output:

enter image description here

.. which is a staircase function but all steps are equal.

Why is my source code not producing the expected output?

Also, I am not sure about my Ntk() function.

Best Answer

Your code already fails within the first function as rexp(0.5) does not generate a realization of an exponential random variable. In fact, the first argument of rexp() is the number of realizations you want to generate (n) followed by the parameter $\lambda$ (rate) as the second argument. So there is no need to define the function rhos() in the first place. In general, your code is quite cumbersome and redundant.

You may simply generate the arrival times as follows

cumsum(rexp(n, lambda)).

or alternatively using the inversion method

cumsum(-log(runif(n))/lambda)

for some specified n and lambda.

Related Question