R/loglikelihood.R

In cbg-ethz/MC-CBN: Large-scale inference of continuous time Conjunctive Bayesian Networks

#' @title Observed Log-Likelihood
#' @export 
#'
#' @description compute the observed log-likelihood
#'
#' @param obs a matrix containing observations or genotypes, where each row
#' correponds to a genotype vector whose entries indicate whether an event has
#' been observed (\code{1}) or not (\code{0})
#' @param poset a matrix containing the cover relations
#' @param lambda a vector of the rate parameters
#' @param eps error rate
#' @param weights an optional vector containing observation weights
#' @param times an optional vector of sampling times per observation
#' @param L number of samples to be drawn from the proposal
#' @param sampling sampling scheme to generate hidden genotypes, \code{X}.
#' OPTIONS: \code{"forward"}, \code{"add-remove"}, \code{"backward"},
#' \code{"bernoulli"}, or \code{"pool"}
#' @param neighborhood.dist an integer value indicating the Hamming distance
#' between the observation and the samples generated by \code{"backward"}
#' sampling. This option is used if \code{sampling} is set to \code{"backward"}.
#' Defaults to \code{1}
#' @param lambda.s rate of the sampling process. Defaults to \code{1.0}
#' @param thrds number of threads for parallel execution
#' @param seed seed for reproducibility
obs.loglikelihood <- function(
  obs, poset, lambda, eps, weights=NULL, times=NULL, L,
  sampling=c('forward', 'add-remove', 'backward', 'bernoulli', 'pool'),
  neighborhood.dist=1L, lambda.s=1.0, thrds=1L, seed=NULL) {
  
  sampling <- match.arg(sampling)
  N <- nrow(obs)
  if (!is.integer(poset))
    poset <- matrix(as.integer(poset), nrow=nrow(poset), ncol=ncol(poset))
  
  if (!is.integer(obs))
    obs <- matrix(as.integer(obs), nrow=N, ncol=ncol(obs))
  
  if (is.null(weights))
    weights <- rep(1, N)

  if (is.null(times)) {
    times <- numeric(N)
    sampling.times.available <- FALSE
  } else {
    lambda.s <- 1 / mean(times)
    sampling.times.available <- TRUE
    if (length(times) != N)
      stop("A vector of length ",  N, " is expected")
  }
  
  if (is.null(seed))
    seed <- sample.int(3e4, 1)
  
  .Call("_obs_log_likelihood", PACKAGE = 'mccbn', obs, poset, lambda, eps,
        weights, times, L, sampling, as.integer(neighborhood.dist), lambda.s,
        sampling.times.available, as.integer(thrds), as.integer(seed))
}

#' @title Complete-data Log-Likelihood
#' @export
#'
#' @description compute the complete-data log-likelihood or, equivalently, the
#' hidden log-likelihood
#'
#' @param lambda a vector of the rate parameters
#' @param eps error rate, \eqn{\epsilon}
#' @param Tdiff a matrix of expected time differences
#' @param dist a vector of expected Hamming distances
#' @param weights an optional vector containing observation weights
complete.loglikelihood <- function(lambda, eps, Tdiff, dist, weights=NULL) {
  
  if (is.null(weights)) {
    weights <- rep(1, length(dist))
  } else if (length(weights) != length(dist)) {
    warning(paste("Length of the weights vector does not agree with the number",
                  " of observations. Setting weights to 1"))
    weights <- rep(1, length(dist))
  }
  
  .Call('_complete_log_likelihood', PACKAGE = 'mccbn', lambda, eps, Tdiff, dist,
        weights)
}