R/logmeanexp.R

In kidusasfaw/pomp: Statistical Inference for Partially Observed Markov Processes

##' The log-mean-exp trick
##'
##' \code{logmeanexp} computes \deqn{\log\frac{1}{N}\sum_{n=1}^N\!e^x_i,}{log
##' mean exp(x_i),} avoiding over- and under-flow in doing so.  It can
##' optionally return an estimate of the standard error in this quantity.
##'
##' When \code{se = TRUE}, \code{logmeanexp} uses a jackknife estimate of the
##' variance in \eqn{log(x)}.
##'
##' @importFrom stats sd
##'
##' @param x numeric
##' @param se logical; give approximate standard error?
##'
##' @return
##' \code{log(mean(exp(x)))} computed so as to avoid over- or
##' underflow.  If \code{se = FALSE}, the approximate standard error is
##' returned as well.
##'
##' @author Aaron A. King
##' 
##' @examples
##'
##' ## an estimate of the log likelihood:
##' ll <- replicate(n=5,logLik(pfilter(ricker(),Np=1000)))
##' logmeanexp(ll)
##' ## with standard error:
##' logmeanexp(ll,se=TRUE)
##'
##' @export
logmeanexp <- function (x, se = FALSE) {
  se <- isTRUE(se)
  mx <- max(x)
  mean <- mx+log(mean(exp(x-mx)))
  if (se) {
    n <- length(x)
    jk <- vapply(seq_len(n),
      function(k) logmeanexp(x[-k]),
      numeric(1))
    c(mean,se=(n-1)*sd(jk)/sqrt(n))
  } else {
    mean
  }
}