# Poisson process simulation I
# STAT 753
# March 12, 2020


############################################################################
# Simulate a homogeneous Poisson process - uses the 'poisson' package in R
############################################################################

# Install poisson package
library("poisson")


####################################################
# Understand the following functions for simulation
####################################################
hpp.event.times
hpp.sim


# Recreate the code above - simulate event times
PP.sim <- function (rate, num.events, num.sims = 1, t0 = 0) 
{
  if (num.sims == 1) {
    x = t0 + cumsum(rexp(n = num.events, rate = rate))
    return(c(t0,x))
  }
  else {
    xtemp = matrix(rexp(n = num.events * num.sims, rate = rate), num.events)
    x = t0 + apply(xtemp, 2, cumsum)
    return(rbind(rep(t0, num.sims), x))
  }
}


# Run the function above, num.sims default is 1
rate = 1
PP.sim(rate, num.events = 10)
PP.sim(rate, num.events = 10, num.sims = 5)

# Compare to hpp.sim
hpp.sim(rate, num.events=10, num.sims = 5)


# Run the function and plot the time series (for 1 simulation)
rate = 1
# rate = 10
num.events = 100
run1 <- PP.sim(rate, num.events)
plot(run1, 0:100, xlab = "Time", ylab = "Number of Events", main = "Poisson Process Sample Path", pch = 19, col="blue")

head(run1)
tail(run1)


#################################################
# Understand the following function for plotting - handles num.sims > 1
#################################################
hpp.plot
?hpp.plot

rate = 1
hpp.plot(rate, num.events = 100)  # default is num.sims=100

hpp.plot(rate, num.events = 100, num.sims = 10)  # default is num.sims=100



###############################################################################
# Simulate a NON-homogeneous Poisson process - uses the 'poisson' package in R
###############################################################################

# Use functions
nhpp.event.times
nhpp.sim

# Similar setup, but requires user to input 'prob.func' which is the non-constant rate (aka intensity function)

# Work through the code for the above 2 functions, run these and plot sample paths for 
# 'prob.func' of your choice

rate <- 1
intensity <- function(t) pmin(t/3, 1)
run2 <- nhpp.sim(rate, num.events = 100, intensity)
plot(run2, 0:100, xlab = "Time", ylab = "Number of Events", main = "Nonhomogeneous Poisson Process Sample Path", 
     pch = 19, col="red")

head(run2)


####################################
# Understand the intensity function
intensity1 = intensity(0:10)
plot(0:10,intensity1, xlab = "Time", ylab = "Intensity (rate)", main = "Intensity Function", 
     type="b", pch = 19, col="darkgreen")

# Beyond time 10, the intensity function is constant at 1
plot(0:20,intensity(0:20), xlab = "Time", ylab = "Intensity (rate)", main = "Intensity Function", 
     type="b", pch = 19, col="darkgreen")


#################################
# Try another intensity function
intensity <- function(t) abs(sin(t))
plot(seq(0,10,by=0.1),intensity(seq(0,10,by=0.1)), xlab = "Time", ylab = "Intensity (rate)", main = "Intensity Function", 
     type="b", pch = 19, col="darkgreen")

# Run the NHPP simulation
rate <- 1
run3 <- nhpp.sim(rate, num.events = 100, intensity)
plot(run3, 0:100, xlab = "Time", ylab = "Number of Events", main = "Nonhomogeneous Poisson Process Sample Path", 
     pch = 19, col="red")

head(run3)



###############################################################################
# Example from R documentation - NHPP is generated by thinning a homogenous PP
###############################################################################

intensity <- function(t) pmin(t/3, 1) # this is a function between 0 and 1
rate <- 10
num.events <- 100
run4 <- nhpp.sim.slow(rate, num.events, prob.func=intensity)
plot(run4, 0:num.events, xlab = "Time", ylab = "Number of Events", 
     main = "Nonhomogeneous Poisson Process Sample Path", pch = 19, col="red")


##################
# Good Example
##################

# Plot intensity function - same example as in NHPP R script for class
lambda <- function(t) (1 + sin(pi*t))
curve(lambda, 0, 10, lwd=2, col="tomato", main="Intensity Function")

# Set constant rate (HPP rate from which to thin)
lambda_max <- 2

# Define prob.func argument as the ratio: lambda/lambda_max
prob.func <- function(t) lambda(t)/lambda_max

# Run the NHPP simulation and plot sample path
rate <- lambda_max
num.events <- 50
run5 <- nhpp.sim(rate, num.events, prob.func=prob.func)
plot(run5, 0:num.events, xlab = "Time", ylab = "Number of Events", 
     main = "Nonhomogeneous Poisson Process Sample Path", pch = 19, col="purple", type="b")