0

I am wondering how to simulate a Poisson process from another one with a Bernoulli random variable of parameter p.

To simulate the first Poisson process with parameter \lambda over the interval [0,t], generally

pois = rpois(1, \lambda)
v = runif(pois, O, t)
w = sort(v)

Now, I know that we can associate a Bernoulli random variable with the arrival times of the first poisson process to simulate another poisson process of parameter p * \lambda but how to do so?

user008
  • 23
  • 4

1 Answers1

0

I think you're asking how to simulate the arrivals time of a Poission process with a given lambda. The time between arrivals in a Poisson process is given by the exponential distribution, so if you want to model x consecutive arrival times of a Poisson process with lamba = 5, you would just do:

cumsum(rexp(x, lambda))

So, for example, suppose I want to model a Poisson arrivals process based on existing Poisson data that I have. (I don't so I'll create a random sample):

set.seed(69)
existing_poisson_data <- rpois(50, 5)
existing_poisson_data
#>  [1]  5  7  6  7  4  8  3  7  3  1  8  4  8  3  3  3  2  6  2  7  4  3  6  4  3  3  6
#> [28]  6  2  4  6  8  4  5  4  6  6  5  4  5 11  4  5  6  3  1  2  3  4  3

I want to simulate the arrivals time for a Poisson process with the same lambda and the same number of arrivals. I can therefore do this:

number_of_arrivals   <- sum(existing_poisson_data)
lambda               <- mean(existing_poisson_data)

simulated_time_diffs <- rexp(number_of_arrivals, lambda)
arrivals             <- cumsum(simulated_time_diffs)

Now we can see whether our arrivals match a Poisson distribution:

simulated_histogram <- hist(arrivals, breaks = 0:ceiling(max(arrivals)))

enter image description here

This looks fine. If the counts are actually Poisson distributed, their mean should be about the same as their variance:

mean(simulated_histogram$counts)
#> [1] 5.418605
var(simulated_histogram$counts)
#> [1] 5.487265

You don't need a dispersion test to see that these data seem to be Poisson distributed. So now you can safely use your arrivals variable as a simulation of the arrival times of a Poisson process with a lambda of 5.

So, if you want to simulate n arrivals for a Poisson process based on another Poisson process with a given lambda and Bernoulli variable p you would do

arrivals <- cumsum(rexp(n, p * lambda))

Alternatively, you could use your Bernoulli variable to sample out some of the arrivals like this:

arrivals <- arrivals[as.logical(rbinom(length(arrivals), 1, p))]
Allan Cameron
  • 147,086
  • 7
  • 49
  • 87
  • thanks for your answer. In fact, there are two ways to simulate the arrivals. The one you did using the exponential, and the second way by ordering the arrivals then assuming that they distribute uniformally so to sample from the uniform. Now my question is: how can I associate the Bernoulli to the arrivals of the first poisson process in order to simulate a newpoisson process of parameter `p * lambda` ? – user008 Feb 03 '20 at 16:33
  • @user008 you can use the way I show on my last line. It assumes that the arrivals will occur with a probability `p` – Allan Cameron Feb 03 '20 at 16:36
  • 1
    @user008 In this case I am simulating a Bernoulli random variable by using rbinom with a size of one. The code then discards any arrivals which don't draw a `1` on this Bernoulli sample. – Allan Cameron Feb 03 '20 at 16:43