Computing log-sum-exponential

log.sum.exp.R

log.sum.exp<- function(x) {
  # Computes log(sum(exp(x))
  # Uses offset trick to avoid numeric overflow: http://jblevins.org/notes/log-sum-exp
  offset <- x[which(x == max(abs(x)))]
  log(sum(exp(x - offset))) + offset
}

And with the possibility to apply only over certain margins/axis:

log.sum.exp <- function(x, margin=NULL) {
  # Computes log(sum(exp(x))
  # Uses offset trick to avoid numeric overflow: http://jblevins.org/notes/log-sum-exp
  innerf <- function(x) {
    if ( max(abs(x)) > max(x) )
      offset <- min(x)
    else
      offset <- max(x)
    log(sum(exp(x - offset))) + offset
  }
  if (is.null(margin)) {
    return(innerf(x))
  }
  else {
    return(apply(x, margin, innerf))
  }
}