archive/likelihood_complete_gsem.R

In sdtemple/incompleteSEM: What the Package Does (One Line, Title Case)

#' Likelihood for complete general stochastic epidemic
#' 
#' Compute exact likelihood for complete epidemic observations.
#' 
#' @param rates numeric vector: rates (beta, gamma)
#' @param r numeric vector: removal times
#' @param i numeric vector: infection times
#' @param N integer: population size
#' 
#' @return negative log likelihood
#' 
#' @export
likelihood_complete_gsem <- function(rates, r, i, N){
  
  beta <- rates[1]
  gamma <- rates[2]
  
  beta <- beta / N
  alpha <- which.min(i)
  n <- length(r)
  
  log_phi <- - beta * (N - n) * sum(r - i)
  log_f <- - gamma * sum(r - i) + n * log(gamma)
  
  psichi <- rep(0, n)
  for(j in (1:n)[-alpha]){
    rj <- r[j]
    ij <- i[j]
    X <- 0
    Y <- 0
    for(k in (1:n)[-j]){
      rk <- r[k]
      ik <- i[k]
      Y <- Y - beta * (min(rk, ij) - min(ik, ij))
      X <- X + as.numeric((ik < ij) & (ij < rk)) * beta
    }
    psichi[j] <- Y + log(X)
  }
  
  log_psichi <- sum(psichi[-alpha])
  #print(-log_phi)
  #print(-log_psichi)
  #print(-log_f)
  
  return(-(log_phi+log_f+log_psichi))
}