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

NHPP_weibull.R

library(animation)

intensity_weibull <- function(t,alpha, beta){
  (beta/alpha)*(t/alpha)^(beta-1)
}

Intensity_weibull <- function(x,alpha,beta){
  (x/alpha)^beta
}

NHPP_weibull = function(Tmax, alpha, beta, maxit){
  mlogz = rexp(1)
  t <- rep(Inf, maxit)
  t[1] = alpha*mlogz^(1/beta[1])
  for(i in 2:maxit){
    mlogz = rexp(1)
    ti = alpha*( mlogz + (t[i-1]/alpha)^beta )^( 1/beta )
    if(ti > Tmax){
      break
    }
    t[i] <- ti
  }
  return(t[1:(i-1)])
}


estimator <- function(dat, t_max){
  n <-length(dat)
  beta <- n/sum(log(t_max)-log(dat))
  alpha <- t_max/(n^(1/beta))
  c("alpha"=alpha, "beta"=beta)
}

###
Tmax <- 50
a <- 1
b <- 1.5
set.seed(1234); ti <- NHPP_weibull(Tmax, a, b, 1000)
estimator(ti, Tmax)

Ts <- seq(5, Tmax, by=1)
saveGIF({
  for( i in seq_along(Ts) ){
    Tmax_c <- Ts[i]
    ti_c <- ti[ti < Tmax_c]
    est <- estimator(ti_c, Tmax_c)
    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")
    curve(Intensity_weibull(x,est[1],est[2]), add=TRUE, col="orange")
    curve(Intensity_weibull(x,a,b), add=TRUE, col="royalblue", lty=2)
  }
}, interval=0.1)