R/mllRH2.R

In RHawkes: Renewal Hawkes Process

mllRH2 <-
function(tms, cens, par,
         h.fn = function(x, p) dexp(x, rate = 1 / p),
         mu.fn = function(x, p){
            exp(dweibull(x, shape = p[1], scale = p[2], log = TRUE) -
            pweibull(x, shape = p[1], scale = p[2], lower.tail = FALSE, log.p = TRUE))
         },
         H.fn = function(x, p) pexp(x, rate = 1 / p),
         Mu.fn=function(x, p){
           - pweibull(x, shape = p[1], scale = p[2], lower.tail = FALSE, log.p = TRUE)
         }){
  p.mu <- par[1:2]
  p.h <- par[3]
  eta <- par[4]
  n <- length(tms)
  if(n == 0) return (Mu.fn(cens, p.mu))
  ## auxiliary functions
  mu <- function(t) mu.fn(t, p.mu)
  Mu <- function(t) Mu.fn(t, p.mu)
  h <- function(t) h.fn(t, p.h)
  H <- function(t) H.fn(t, p.h)
  phi <- function(s) eta * sum(h(s - tms[tms < s]))
  Phi <- function(t) eta * sum(H(t - tms[tms < t]))
  ## vectors to hold log(p(tau_i|tau_{1:i-1})), log(d_{i,1:(i-1)}),
  ## log(p_{i,1:(i-1)}), log(p_{i-1,1:(i-2)})
  lden <- lcon.den <- numeric(n);
  lpi <- lpimo <- numeric(n);
  A <- U <- numeric(n); 
  A[1] <- Mu(tms[1]); U[1] <- 1-exp(-A[1])
  res <- 0;
  lden[1] <- -Mu(tms[1]) + log(mu(tms[1]));
  res <- res - lden[1]
  i <- 2;
  while(i <= n + 1){
    if(i == 2){
      lpi[1] <- 0;
      if(i <= n){
        lcon.den[1] <- - Mu(tms[2] - tms[1]) - Phi(tms[2]) +
          log(mu(tms[2] - tms[1]) + phi(tms[2]))
      }else{
        lcon.den[1] <- - Mu(cens - tms[1]) - Phi(cens)
      }
    }else{
      ph <- phi(tms[i - 1])
      m <- mu(tms[i - 1] - tms[1:(i - 2)])
      lpi[1:(i - 2)] <-
        log(ph) - log(ph + m) + lcon.den[1:(i - 2)] + 
        lpimo[1:(i - 2)]- lden[i - 1]
      mx <- max(lpimo[1:(i - 2)] + lcon.den[1:(i - 2)])
      lpi[i - 1] <- log(weighted.mean(m / (ph + m),
                        exp(lcon.den[1:(i - 2)] + lpimo[1:(i - 2)] - mx))
      )
      if(i <= n){
        lcon.den[1:(i - 1)] <- - Mu(tms[i] - tms[1:(i - 1)]) +
          Mu(tms[i - 1] - tms[1:(i - 1)]) -
          Phi(tms[i]) + Phi(tms[i - 1]) +
          log(mu(tms[i] - tms[1:(i - 1)]) + phi(tms[i]))
      }else{
        lcon.den[1:(i - 1)] <- - Mu(cens - tms[1:(i - 1)]) +
          Mu(tms[i - 1] - tms[1:(i - 1)]) -
          Phi(cens) + Phi(tms[i - 1])
      }
    }
    mlcon.den <- max(lcon.den[1:(i - 1)])
    mlpi <- max(lpi[1:(i - 1)])
    lden[i] <- log(weighted.mean(exp(lcon.den[1:(i - 1)] - mlcon.den),
                    exp(lpi[1:(i - 1)] - mlpi))) + mlcon.den
    res <- res - lden[i]
    lpimo[1:(i - 1)] <- lpi[1:(i - 1)]
    tmp <- eta * (sum(H(tms[i] - tms[tms < tms[i]])) - 
              sum(H(tms[i - 1] - tms[tms < tms[i - 1]])))
    A[1:(i - 1)] <- Mu(tms[i] - tms[tms < tms[i]]) - 
                        Mu(tms[i - 1] - tms[tms < tms[i]]) + tmp
    U[i] <- sum(exp(lpi[1:(i - 1)]) * (1 - exp(- A[1:(i - 1)])))
    i <- i + 1
  }
  U <- U[-(n + 1)] #there is no (n + 1)th residual
  list(mll = res,U = U,n = n)
}