# Nonhomogeneous Poisson process simulation 
# STAT 753
# March 24, 2020


#############################################
# Simulate a non-homogeneous Poisson process
#############################################

NHPP.sim <- function(lambda, t_max=10){
  #set.seed(1)
  t <- 0
  # Constant rate that is the max rate given by the intensity function lambda
  lambda_star <- function() max(sapply(seq(1, t_max,length.out=1000), lambda))
  X <- numeric()
  while(t <= t_max){
    u <- runif(1)
    # Update time point
    t <- t - log(u)/lambda_star()
    # Probability of acceptance of the new time point
    if(runif(1) < lambda(t)/lambda_star()) {
      X <- c(X,t)
    }
  }
  return(c(0,X))
}


###########################
# Example 1
###########################

# Intensity function: 
#   lambda(t) = 1 + b*t for different values of b (see below)

# Plot intensity function for b=1, for example
b <- 1
lambda <- function(t) 1 + b*t
curve(lambda, 0, 10, lwd=2, main="Intensity Function")


# Set max time to run the simulation
t_max = 20
  
# Run 1
b <- 1
lambda <- function(t) 1 + b*t
run1 <- NHPP.sim(lambda,t_max)
n1 <- length(run1)-1

# Run 2
b <- 0.1
lambda <- function(t) 1 + b*t
run2 <- NHPP.sim(lambda,t_max)
n2 <- length(run2)-1

# Run 3
b <- 0.01
lambda <- function(t) 1 + b*t
run3 <- NHPP.sim(lambda,t_max)
n3 <- length(run3)-1

# Find longest sample path and max height
nmax <- max(n1,n2,n3)
hmax <- ceiling(max(run1,run2,run3))

# Plot the 3 sample paths on the same graph
plot(run1, 0:n1, type="b", pch=19, col="blue", xlim=c(0,hmax), ylim=c(0,nmax), xlab="Time", ylab="N(t)",
     main="Nonhomogeneous Poisson Process: lambda(t) = 1+bt") 
points(run2, 0:n2, type="b", pch=19, col="red")
points(run3, 0:n3, type="b", pch=19, col="aquamarine4")
legend("topleft", c("b=1","b=0.1","b=0.01"), pch=c(19,19,19), col=c("blue","red","aquamarine4"))



###########################
# Example 2
###########################

# Intensity function: 
#   lambda(t) = (1 + sin(pi*t))

# Plot intensity function
lambda <- function(t) (1 + sin(pi*t))
curve(lambda, 0, 10, lwd=2, col="tomato", main="Intensity Function")


# Set max time to run the simulation
t_max = 50

# Run 1
run1 <- NHPP.sim(lambda,t_max)
n1 <- length(run1)-1

# Run 2
run2 <- NHPP.sim(lambda,t_max)
n2 <- length(run2)-1

# Run 3
run3 <- NHPP.sim(lambda,t_max)
n3 <- length(run3)-1

# Find longest sample path and max height
nmax <- max(n1,n2,n3)
hmax <- ceiling(max(run1,run2,run3))


# Plot the 3 sample paths on the same graph
plot(run1, 0:n1, type="b", pch=19, col="blue", xlim=c(0,hmax), ylim=c(0,nmax), xlab="Time", ylab="N(t)",
     main="Nonhomogeneous Poisson Process: lambda(t) = 1+sin(pi*t)") 
points(run2, 0:n2, type="b", pch=19, col="red")
points(run3, 0:n3, type="b", pch=19, col="aquamarine4")
legend("topleft", c("run 1","run 2","run 3"), pch=c(19,19,19), col=c("blue","red","aquamarine4"))

# Save one run to plot with Example 3 below
#lambda <- function(t) (1 + sin(pi*t))
run1 <- NHPP.sim(lambda,t_max=10)
run_example2 <- run1



###########################
# Example 3
###########################

# Modify the intensity function from Example 2 - e.g. increase/decrease the frequency and amplitude of lambda
# and see how it affect the NHPP sample paths, compare to version given in Example 2.

# Intensity function: 
#   lambda(t) = a*(1 + sin(b*pi*t))

# Plot intensity function (a>1 increases the amplitude, b>1 increases the frequency)
a <- 5  # or 10, 20
b <- 1   # or 2

lambda <- function(t) 1 + a*(1 + sin(b*pi*t))
curve(lambda, 0, 10, lwd=2, col="purple", main="Intensity Function")


# Set max time to run the simulation
t_max = 10

# Run 1
run1 <- NHPP.sim(lambda,t_max)
n1 <- length(run1)-1

# Run 2
run2 <- NHPP.sim(lambda,t_max)
n2 <- length(run2)-1

# Run 3
run3 <- NHPP.sim(lambda,t_max)
n3 <- length(run3)-1

# Find longest sample path and max height
nmax <- max(n1,n2,n3)
hmax <- ceiling(max(run1,run2,run3))


# Plot the 3 sample paths on the same graph
plot(run1, 0:n1, type="b", pch=19, col="blue", xlim=c(0,hmax), ylim=c(0,nmax), xlab="Time", ylab="N(t)",
     main="Nonhomogeneous Poisson Process: lambda(t) = 1+a(1+sin(b*pi*t))") 
points(run2, 0:n2, type="b", pch=19, col="red")
points(run3, 0:n3, type="b", pch=19, col="aquamarine4")
#legend("topleft", c("b=1","b=0.1","b=0.01"), pch=c(19,19,19), col=c("blue","red","aquamarine4"))

# Compare to run from Example 2 - if t_max differs, that's ok
points(run_example2, 0:(length(run_example2)-1), type="b", pch=19, col="orange")
legend("topleft", c("run 1","run 2","run 3", "run ex2"), pch=c(19,19,19,19), 
       col=c("blue","red","aquamarine4","orange"))