# Nonhomogeneous poisson process simulation

I've been looking at ways to generate a Nonhomogeneous Poisson Process (NHPP) including the nonlinear time transformation (using a rate-1 process and inverting the cumulative rate function). I've also seen the thinning approach.

QUESTION: How do you generate the NHPP via thinning? What are the key differences between the two approaches? What would make one prefer one method over the other?

CONTEXT I am trying to use a non-homogeneous poisson process to simulate the claims next year, the rate varies depending on which month, my rate function is:

nhpp_lambda <- function(t) {
if (t <= 59) {lambda = 20.83}
else if (t > 59 & t <= 151) {11.02}
else if (t > 151 & t <= 243) {11.68}
else if (t > 243 & t <= 334) {26.41}
else if (t > 334 & t <= 365) {20.83}
}


When I try to adopt the thinning approach and run the following code:

nhpp_simulate <- function() {
d <- 0
X <- 0
while(X[length(X)] <= 365) {
U <- runif(1)
d <- d - log(U)/lambda_star
U2 <- runif(1)
if(U2 <= nhpp_lambda(d)/27) {X <- c(X,d)}
}
return(X)
}

blahblah <- nhpp_simulate()


R gives the error "Error in if (U2 <= nhpp_lambda(d)/27) { : argument is of length zero" but length(U2) is definitely not zero?

• which packages does nhpp_lambda belong to please? It is not NHPP is it? Commented Jan 23, 2022 at 14:28

Thinning Approach: (method below by OP's request)
Inefficient when fluctuation in time is large (big $$\bar \lambda$$ gives a high rejection probability). Possible workaround is to break interval $$[0,T]$$ into small intervals and pick a $$\bar \lambda$$ for each interval.

Nonlinear Time Transformation (inversion) Approach: see method here
Inefficient when $$m(t)$$ is difficult to invert (requiring numerical search). Possible workaround is to consider a piecewise constant approximation for $$\lambda(t)$$.

Generating a NHPP: Thinning Method
Task: Generate a NHPP $$\{N(t),t\ge 0\}$$ with rate $$\lambda(t)$$ on the interval $$t\in [0,T]$$.
Preparation:
Choose $$\bar \lambda$$ such that $$\lambda(t) \le \bar \lambda$$ for all $$t\in [0,T]$$.
Define the thinning probability $$p(t) = \frac{\lambda(t)}{\bar\lambda}$$.

Procedure:
Let $$I$$ be the number of events that occur by time $$t$$. Then $$S(I)$$ is the time of the most recent event.
1. $$t=0$$;$$\,I = 0$$.
2. Generate $$U_1\sim \text{Uniform}(0,1)$$.
3. Set $$t = t - (\text{ln}(U_1)/\bar \lambda)$$. If $$t>T$$, stop; else go to 4.
4. Generate $$U_2 \sim \text{Uniform}(0,1)$$, independent of $$U_1$$.
5. If $$U_2 \le p(t)$$, set $$I = I + 1$$; $$S(I) = t$$.
6. Go to step 2.

At the end of the procedure, you have event times and the counting process.

Validation Since for a fixed $$t$$, $$N(t)\sim\text{Poisson}(m(t))$$ with mean $$m(t) = \int_0^t \lambda(s)ds$$, it is easy to validate the results.
1. Ensure the Dispersion equals 1.
2. Ensure you match the rate function (arrival or cumulative).

Let $$IDC(t) = \frac{\text{Var}(N(t))}{\text{E}[N(t)]}$$ which is the Index of Dispersion (for counts). For any NHPP, $$IDC(t) = 1,\, \forall t$$.

The figure below demonstrates the correctness of the procedure.

The target arrival rate and cumulative arrival rate are given below.

Bootleg MATLAB code...

function [ EventTimes ] = genNHPP(rate_fh,T,n,dt)
%GENNHPP Generate a NHPP
%   Generates a NHPP, N(t): # of events by time t
% INPUTS
%   rate_fh : function handle for rate (vectorized)
%   T : end of time horizon
%   n : number of sample paths to generate (default n = 1)
% OUTPUTS
%   EventTimes : if n = 1, EventTimes is a vector (length is number of events)
%                if n > 1, EventTimes is a cell array with n rows (# columns is number of events)
%   NumArrived : Number Arrived by time t (time is row, columns are sample paths)

% Input Error Checking ****************************************************
narginchk(2, 4)
if nargin < 3 || isempty(n), n = 1; end % Default is 1 sample path
if nargin < 4 || isempty(dt), dt = 1; end % Default is 1 day
if ~isa(rate_fh,'function_handle')
error('rate_fh must be a valid function handle')
end
if T <= 0, error('T must be a positive real number'), end
% End (Input Error Checking) **********************************************

% Generate a NHPP
MaxLambda = max(rate_fh([0:10^-5:T]));
ph=@(t) rate_fh(t)/MaxLambda;
if n ==1
% Single Sample Path
t = 0; Nevents = 0; EventTimes = []; done = false;
while ~done
t = t + (-1/MaxLambda)*log(rand());
if t > T
done = true;
else
if rand()<= ph(t);
Nevents = Nevents+1;
EventTimes(Nevents) = t;
end
end
end

else
% Multiple Sample Paths
NumPaths = n; EventTimes = {};
for path = 1:NumPaths
t = 0; Nevents = 0; done = false;
while ~done
t = t + (-1/MaxLambda)*log(rand());
if t > T
done = true;
else
if rand() <= ph(t);
Nevents = Nevents+1;
EventTimes{path,Nevents} = t;
end
end
end
end
end

• Thank you for the detailed answer!! I am not extremely familiar with coding in matlab, but i read through the steps in the procedure you provided, i think thats what i tried to do myself, is something wrong with my code? Commented Sep 29, 2018 at 6:31
• if I change my rate function to: nhpp_lambda <- function(t) { if (t > 0 && t <= 59) {lambda = 20.83} else if (t > 59 && t <= 151) {lambda = 11.02} else if (t > 151 && t <= 243) {lambda = 11.68} else if (t > 243 && t <= 334) {lambda = 26.41} else if (t > 334 && t <= 365) {lambda = 20.83} else if (t > 365) {lambda = 0} return(lambda) } then the code can run, am I allowed to do this? Commented Sep 29, 2018 at 8:50
• I'm unsure what you mean by are you 'allowed'? What is length of nhpp_lambda(d)/27 ? I'm not very good with R so I'm not sure how to help. Maybe try a script instead of a function and manually increment loop so you can evaluate step by step and see where it goes wrong. Commented Sep 29, 2018 at 14:24
• If i dont add that extra bit to my rate function, for some reason my d becomes larger than 365 Commented Sep 29, 2018 at 21:42