Anime for the point estimation of non-homogenous Poisson process with exponentially growth

NHPP_exp.R

library(animation)

#累積強度関数
Intensity_exp <- function(t, a, b) {
  ##  (exp(a*t+b)-exp(b))/a
  exp(b) * expm1(a * t) / a
}

#強度関数
intensity_exp <- function(t,a,b){
  exp(a*t+b)
}

loglik_exp <- function(par, t, Tmax){
  a <- par[1]
  b <- par[2]
  - sum(a*t+b) + Intensity_exp(Tmax,a,b)
}

grad_loglik_exp <- function(par, t, Tmax){
  a <- par[1]
  b <- par[2]
  ga <- - sum(t) + exp(b)*( (exp(a*Tmax)*(a*Tmax-1)+1)/a^2 )
  gb <- - length(t) + Intensity_exp(Tmax,a,b)
  c(ga, gb)
}

NHPP_exp = function(Tmax, a, b, maxit){
  mlogz = rexp(1)
  t <- rep(Inf, maxit)
  t[1] <- (log(a*mlogz + exp(b))-b)/a
  for(i in 2:maxit){
    mlogz = rexp(1)
    s <- a*mlogz + exp(a*t[i-1]+b) 
    if(s<0){
      break
    }
    ti =(log(s) - b)/a
    if(ti>Tmax){
      break
    }
    t[i] <- ti
  }
  return(t[1:(i-1)])
}


#####
Tmax <- 100
a <- 0.02
b <- 0
set.seed(1234); ti <- NHPP_exp(Tmax, a, b, 1000)
Ts <- seq(5, 100, by=2)
saveGIF({
  for( i in seq_along(Ts) ){
    Tmax_c <- Ts[i]
    ti_c <- ti[ti < Tmax_c]
    opt1 <- optim(c(a,b), loglik_exp, gr = grad_loglik_exp,
                  t=ti_c, Tmax=Tmax_c, 
                  method = "BFGS", hessian = TRUE)
    
    plot(c(0,ti), c(0,seq_along(ti)), type="s", xlab="time", ylab = "cumulative count", col="grey")
    lines(c(0,ti_c), c(0,seq_along(ti_c)), type="s")
    abline(v=Tmax_c, lty=3)
    curve(Intensity_exp(x,opt1$par[1],opt1$par[2]), add=TRUE, col="orange")
    curve(Intensity_exp(x,a,b), add=TRUE, col="royalblue", lty=2)
  }
}, interval=0.1)