# How to simulate Poisson arrival times if the rate varies with time?

Suppose we are to study a non-homogeneous Poisson process of 3 hour cycles in which:

At the first hour, the arrival rate is 1.5 events / hr.

At the second hour, the arrival rate is 2.1 events / hr.

At the third hour, the arrival rate is 3.4 events / hr.

I would like to generate the time points (not inter-arrival times) of events observed in $3\times1000$ hours. How should I implement this in R?

Edit1: Following jbowman's answer, I tried to verify the mean rate with

ia=diff(arrival_times)
t=arrival_times[2:length(arrival_times)]%%3
ia_bygrp=list(ia[which(0<=t & t<1)],ia[which(1<=t & t<2)],ia[which(2<=t & t<3)])
sapply(ia_bygrp,function(x) 1/mean(x))


The result was, unexpectedly (I only quote 1 result but I tried for 30+ times to reduce randomness)

[1] 1.892691 2.869797 2.116035


Why? See Edit2.

Edit2: I tried to use larger rates: rates = c (15.0, 21.0, 34.0) for example. The arithmetic means were instead

[1] 21.64794 32.72656 15.59476


I believe it's because of the small number of occurrence n_arrivals in each hour which makes runif(n_arrivals) not-so-uniform. Any other solution?

Edit3: Based on GeoMatt22's code with Octave (well, there is large variety here), I wrote another code to simulate the empirical rates $\lambda$.

Here is the code

rates <- c(1.5, 2.1, 3.4)
arrival_times <- c()
Nrep=10000
for (hour in 1:(3*Nrep)) {
n_arrivals <- rpois(1,rates[1 + hour%%3])
arrival_times <- c(arrival_times, runif(n_arrivals) + hour - 1)
}
arrival_times <- sort(arrival_times)

x=hist(arrival_times%%3,breaks=c(0,1,2,3),plot=FALSE)
bp=barplot(x$counts/Nrep,axes=FALSE,ylim=range(0,max(rates)*1.1),ylab="Observed rate / hr",main="Simulation for jbowman's answer") points(bp,rates[c(2,3,1)],col='red',pch=20) axis(side=2,at=seq(0,3.5,0.1),label=seq(0,3.5,0.1)) axis(side=1,at=bp,label=c("1st hours","2nd hours","3rd hours"))  and the histogram (red dot = true rates) Pretty close. Apparently counting occurrences in each hour is not the same as the reciprocal of mean hourly inter-arrival times. [I think I have seen a question about the difference between the two but I could not retrieve it at the moment]. It turns out jbowman's answer can simulate well even using rates down to rates <- c(1.5, 2.1, 3.4)/100. • This is a coding question. I don't think OP has an issue with the conceptual understanding of the problem at all – Aksakal Apr 29 '17 at 13:52 • @Aksakal the main question is about simulating a piecewise-homogeneous Poisson process. There are some R aspects in the edits, but I think the core question is statistical. (IMO this is worth keeping, if it is not a duplicate.) – GeoMatt22 Apr 29 '17 at 14:13 • @GeoMatt22, you are right, I'm retracting my close vote – Aksakal Apr 29 '17 at 15:22 ## 5 Answers There are already several reasonable answers, but to try and add some clarity, here I will try to synthesize these, and do some verification on the resulting algorithm. Homogeneous Process For a homogeneous 1D Poisson process with constant rate$\lambda$, the number of events$n$in a time interval$T$will follow a Poisson distribution $$n\sim\text{Poiss}_\lambda$$ And the arrival times$t_1,\ldots,t_n$will be i.i.d. uniformly through the interval $$t_i\sim\text{Unif}_{[0,T]}$$ while the waiting times between events will be i.i.d. exponentially $$\tau_i=t_{i+1}-t_i\sim\text{Exp}_\lambda$$ with CDF $$\Pr\big[\tau<\Delta{t}\big]=1-\exp^{-\lambda\Delta{t}}$$ So a homogeneous Poisson process can be easily simulated by first sampling$n$and then sampling$t_{1:n}$(or, alternatively sampling$\tau$until$t=\sum\tau>T$). Inhomogeneous Process For an inhomogeneous Poisson process with rate parameter$\lambda(t)$the above can be generalized by working in the transformed domain $$\Lambda(t)=\int_0^t\lambda(s)ds$$ where$\Lambda(t)$is the expected cumulative number of events up to time$t$. (As noted in Aksakal's answer, and the reference cited in jth's answer.) In the generalized approach we simulate in the$\Lambda$space, where "time" is dilated so in the deformed timeline the Poisson process is now a homogeneous process with unit rate. After simulating this$\Lambda$process, we map the samples$\Lambda_{1:n}$into arrival times$t_{1:n}$by inverting the (monotonic)$\Lambda(t)$mapping. Note that for the piecewise-constant rates$\lambda(t)$here, the mapping$\Lambda(t)$is piecewise linear, so very easy to invert. Example Simulation I made a short code in Octave to demonstrate this (listed at the end of this answer). To try and clear up questions about the validity of the approach for small rates, I simulate an ensemble of concatenated simulations. That is, we have 3 rates$\boldsymbol{\lambda}=[\lambda_1,\lambda_2,\lambda_3]$each of duration 1 hour. To gather better statistics, I instead simulate a process with$\boldsymbol{\hat{\lambda}}=[\boldsymbol{\lambda},\boldsymbol{\lambda},\ldots,\boldsymbol{\lambda}]$, repeatedly cycling through the$\boldsymbol{\lambda}$vector to allow larger sample sizes for each of the$\lambda$'s (while preserving their sequence and relative durations). The first step produces a distribution of waiting times (in$\Lambda$space) like this which compares very well to the expected theoretical CDF (i.e. a unit-rate exponential). The second step produces the final ($t$space) arrival times. The empirical distribution of waiting times matches well with theory here as well: This figure is a little bit more busy, but can be interpreted as follows. First, the upper panel shows the component CDFs associated with a homogeneous Poisson process at each of the three$\lambda$'s. We expect the aggregate CDF to be a mixture of these components $$\Pr\big[\tau<\Delta{t}\big]=\sum_iw_i\Pr\big[\tau<\Delta{t}\mid\lambda=\lambda_i\big]$$ where$w_i=\Pr\big[\lambda=\lambda_i\big]$are the mixing fractions. For each component process, the expected number of samples is$\langle{n_i}\rangle=\lambda_iT_i$. Since the durations$T_i$are equal, the expected mixing weights will then scale with the rates, i.e. $$w_i\propto\lambda_i$$ The lower panel above shows the expected CDF obtained by mixing the components in proportion to their rates. As can be seen, the empirical CDF of the inhomogeneous simulation is consistent with this theoretically expected CDF. Example Code The following code can be run here. (Note: the algorithm is only 4 lines; the bulk of the code is comments and verification.) %% SETUP lam0=[1.5,2.1,3.4]; dt0=ones(size(lam0)); % rates and their durations Nrep=1e3; lam=repmat(lam0,[1,Nrep]); dt=repmat(dt0,[1,Nrep]); % do replications (for stats check) %% SIMULATION L=cumsum([0,lam.*dt]); t=cumsum([0,dt]); % cumulative expected # events and time Lmax=L(end); N=poissrnd(Lmax); % sample total # events Lsmp=Lmax*sort(rand([1,N])); % sample event times (in "L space") tsmp=interp1(L,t,Lsmp); % transform to "t space" %% STATS CHECK % "L space" waiting time CDF dL=sort(diff(Lsmp)); p=(1:N-1)/N; % simulated p0_L=1-exp(-dL); % exponential h=plot(dL,p,'k',dL,p0_L,'r','LineWidth',1.5); set(h(1),'LineWidth',4); axis tight; xlabel('dL'); legend('data','theory (exp)',4); title('L space CDF'); % "t space" waiting time CDF dT=sort(diff(tsmp)); % simulated wcmp=(lam0.*dt0)/(lam0*dt0'); pcmp=1-exp(-lam0'*dT); p0=wcmp*pcmp; % mixture subplot(211); plot(dT,pcmp); ylabel('CDF'); title('Mixture Components'); axis tight; legend(cellstr(num2str(lam0','lam = %0.1f ')),4); subplot(212); h=plot(dT,p,'k',dT,p0,'r','LineWidth',1.5); set(h(1),'LineWidth',4); axis tight; xlabel('dt'); ylabel('CDF'); title('Aggregate'); legend('data','theory (exp mixture)',4); % mean arrival rate tbin=0.5:2.5; Navg=hist(mod(tsmp,3),tbin)/Nrep; bar(tbin,Navg,1); hold on; plot(tbin,lam0,'r.'); hold off; legend('observed (hist.)','theory (lam.)',0); xlabel('hour'); ylabel('arrivals/hour');  As requested in the comments, here is an example of the simulated mean arrival rates (see also updated code above): • I wonder if we could include calculate the empirical rate$\lambda_i$from$\frac{1}{dT_i}$(segregated into different time sessions)? – user2513881 Apr 30 '17 at 12:03 • You just do a histogram on arrival times modulo 3 (see update). If the answer was helpful, feel free to upvote it :) BTW the idea of concatenating the simulations (to get larger samples) should apply to any of the techniques suggested in the other answers as well. – GeoMatt22 Apr 30 '17 at 15:41 • I decided to accept this answer as it illustrated fully the conceptual and statistical aspects (even though the answer originated from others). The final histogram just made the whole answer complete to me. – user2513881 Apr 30 '17 at 15:52 It's pretty straightforward, thanks to the fact that the interarrival times of a Poisson process are exponentially distributed and the exponential distribution has the memoryless property. (As an aside, note that if you are generating the time points, you are generating the interarrival times too, as they are just the differences between successive time points, and vice versa.) One strategy is to generate the number of arrivals in each of the 3 x 1000 hours, then make use of the fact that the arrival times are uniformly distributed within each hour (thanks to the memoryless property.) If you have, say, 3 arrivals in hour 1, you can generate the arrival times for that hour by generating three uniform variates over the hour. You can ignore "crossing the hour boundary" effects, also because of the memoryless property of the exponential distribution. For example: rates <- c(1.5, 2.1, 3.4) arrival_times <- c() for (hour in 1:3000) { n_arrivals <- rpois(1,rates[1 + hour%%3]) arrival_times <- c(arrival_times, runif(n_arrivals) + hour - 1) } arrival_times <- sort(arrival_times)  The first and last few arrival times are: > c(head(arrival_times), tail(arrival_times)) [1] 6.175204e-02 5.907350e-01 1.066275e+00 1.089332e+00 1.638492e+00 2.296899e+00 2.996753e+03 2.996817e+03 2.997005e+03 [10] 2.997376e+03 2.998619e+03 2.999689e+03  • I tried to verify the mean of interarrival times. But the result is not as expected. Please see the edited question. – user2513881 Apr 29 '17 at 4:34 • @user2513881 for such small numbers of samples, it is expected that the mean will vary a lot. To judge you could check the associated confidence intervals (e.g. see here). Or, more simply, just simulate a bunch of times (e.g. until$\sum n_\text{event}$exceeds some threshold) and look at the ensemble statistics. – GeoMatt22 Apr 29 '17 at 5:42 • @GeoMatt22: "I only quote 1 result but I tried for 30+ times to reduce randomness". I further tried 10,000 simulations (each with 3,000 hours). The aggregated mean turns out to be [1] 2.473990 2.587038 1.797931. For the hours simulated with rate 3.4, the range is 2.268 - 2.906. Obviously CI will not cross the expected mean of 3.4. Running the simulation with 1 arrival for 1000 times is not the same as running another with 100 arrival for 10 times. This simulation algorithm does not converge to mean rate for small rates, seemingly. – user2513881 Apr 29 '17 at 6:07 • @user2513881 OK, that makes sense. But does that mean the simulation algorithm is wrong? For a low rate and a short timespan, the exponential will necessarily have a severely truncated tail. Try the following: First simulate 1000 times, i.e.$\vec{t}_i$for$i=1:1000$. Then concatenate the event times together, offset by$i-1$, i.e.$[\vec{t}_1,\vec{t}_2+1,\ldots ]$. When I do this, the concatenated$\Delta\vec{t}$has the appropriate exponential PDF. – GeoMatt22 Apr 29 '17 at 7:43 • In Octave I did this: lam=3.4; N=1e3; T=[]; for t0=1:N; n=poissrnd(lam); t=sort(rand(1,n)); T=[T,t+t0-1]; end; dT=diff(T); m=numel(dT); p=(1:m)/(m+1); dT=sort(dT); plot(dT,p,dT,1-exp(-lam*dT)); – GeoMatt22 Apr 29 '17 at 7:44 There are several methods. For simple rate functions like yours I find the inversion approach to be the easiest. See algorithm 3 of this paper; it's 5 lines. Edit: and the algorithm is: rates = c(1,2,3) s = sum(rates) cs = c(0, cumsum(rates)) R = function(t) { L = length(rates) t0 = floor(t) d = t0 %% L + 1 return(s * floor(t / L) + cs[d] + (t - t0) * rates[d]) } Rinv = function(y) { t0 = L * y / s return(uniroot(function(t) R(t) - y, c(t0, t0 + L))$root)
}
# Actual algorithm itself is only 2 lines.
x = cumsum(rexp(1000))  # desired nbr of points
y = sapply(x, Rinv)


If efficiency matters then you'd want to analytically invert R(t).

• Can you add the 5 lines? (To make your answer self-contained and avoid link-rot, see here.) – GeoMatt22 Apr 29 '17 at 7:51
• – GeoMatt22 Apr 29 '17 at 8:03

Your intensities are so low that it's impossible to ignore the boundary effects. By that I mean that a time between the last event in the first hour and the first event in the next hour is covered by two Poisson (exponential) processes.

Method 1: So, the simplest way to treat this (certainly, not the most efficient) is to construct the nonhomegenous intensity $\lambda(t)$ function, where $t$ is time. Next, step in time by one second (or a smaller interval), then generate the Bernoulli experiments, until you get to the end of the time.

Method 2: However, a much better (faster) solution is to scale the time instead of scaling the intensity. Here's the idea. You have a homogenous Poisson process, BUT the time speeds up after one hour, NOT the intensity increases. Your clock ticks faster after the first hour by a factor of $\frac{2.1}{1.5}=1.4$ or $\frac{3.4}{1.5}\approx 2.27$. So, all you need to do is to generate a bunch of exponential arrival times, then compress the times after the first hour by appropriate factors.

Solution for Method 1

This is not R forum, and other folks might be interested in the solution, so here's a pseudo code:

function y = lambda(t)
% nonhomegenous Poisson intensity
% t - time is hours
if t<1
y = 1.5;
else if t<2
y = 2.1;
else
y = 3.4;
end if
end

function times = f()
% times in hours of the event occurrences
list times;
dt = 1/3600;
for t = 0 to 3 step dt
if rand() < lambda(t)*dt
end if
end
end


Solution for Method 2

function times = f()
lambda = 1.5
list times
t = 0
do
dt = random_exponential(lambda)
t = t + dt
st = compress(t)
if st<3
end if
while st < 3
end

% compress slow time to faster times after 1st hour
function st = compress(t)
scale2 = 2.1/1.5
scale3 = 3.4/1.5
t2 = 1 + scale2;
st = min(t,1) + (min(t,t2)-min(t,1))/scale2 + (max(t,t2)-t2)/scale3
end

• (+1) I think this is the same approach as @jth, but the added explanation and generic pseudocode is much more clear. – GeoMatt22 Apr 29 '17 at 15:25
• If dt > 3 in Method 2 then is the calculated rate correct? – jth Apr 29 '17 at 15:43
• @jth, that's a bug, you're right – Aksakal Apr 29 '17 at 15:54
• I bit the bullet and tried my own hand at an answer, based on the time-dilation idea suggested by you and @jth. About the same, but I tried to add some verification to hopefully increase the OP's confidence in the result (and increase the non-R code diversity, I guess :) ). – GeoMatt22 Apr 29 '17 at 18:54
• You are trying to do this by generating interarrival times, but the OP said he/she did not want to generate interarrival times. Also, you can ignore the boundary effects if you know how many events occurred in each hour, and you can find that out by randomly generating them, because they are distributed Poisson with the hour-specific rate. – jbowman Apr 29 '17 at 21:06

If Python is not a problem for you, I advise you to have a look at this new open-source library: tick. Among others, it provides simulation and inference methods for point processes. You can find an example explaining the simulation of an inhomogeneous Poisson process here.

• Thanks for your suggestion but I am more inclined to use R. – user2513881 Apr 28 '17 at 13:49