library("matrixStats")
library("stats")
logSumExp_R <- function(lx, na.rm = FALSE) {
log(sum(exp(lx), na.rm = na.rm))
}
## R-help thread \emph{'[R] Beyond double-precision?'} on May 9, 2009.
for (mode in c("integer", "double")) {
cat("mode: ", mode, "\n", sep = "")
set.seed(1)
x <- runif(20, min = 1.0, max = 3.0)
storage.mode(x) <- mode
str(x)
## The logarithm of the harmonic mean
y0 <- log(1 / mean(1 / x))
print(y0) ## -1.600885
lx <- log(x)
y1 <- log(length(x)) - logSumExp(-lx)
print(y1) ## [1] -1.600885
# Sanity check
stopifnot(all.equal(y1, y0))
y2 <- log(length(x)) - logSumExp_R(-lx)
# Sanity check
stopifnot(all.equal(y2, y0))
} # for (mode ...)
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# Missing values
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
## NA values
lx <- c(1, 2, 3)
lx[2] <- NA_real_
y0 <- logSumExp_R(lx, na.rm = FALSE)
y <- logSumExp(lx, na.rm = FALSE)
print(y)
stopifnot(identical(y, NA_real_))
stopifnot(all.equal(y, y0))
y0 <- logSumExp_R(lx, na.rm = TRUE)
y <- logSumExp(lx, na.rm = TRUE)
print(y)
stopifnot(all.equal(y, y0))
## NaN values
lx <- c(1, 2, 3)
lx[2] <- NaN
y0 <- logSumExp_R(lx, na.rm = FALSE)
y <- logSumExp(lx, na.rm = FALSE)
print(y)
stopifnot(identical(y, NA_real_))
stopifnot(all.equal(y, y0))
y0 <- logSumExp_R(lx, na.rm = TRUE)
y <- logSumExp(lx, na.rm = TRUE)
print(y)
stopifnot(all.equal(y, y0))
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# Corner cases
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
## Zero-length vectors
lx <- numeric(0L)
y0 <- logSumExp_R(lx)
y <- logSumExp(lx)
print(y)
stopifnot(identical(y, -Inf))
stopifnot(all.equal(y, y0))
## Vectors of length one
lx <- 1.0
y0 <- logSumExp_R(lx)
y <- logSumExp(lx)
print(y)
stopifnot(identical(y, lx))
stopifnot(all.equal(y, y0))
lx <- NA_real_
y0 <- logSumExp_R(lx, na.rm = TRUE)
y <- logSumExp(lx, na.rm = TRUE)
print(y)
stopifnot(identical(y, -Inf))
stopifnot(all.equal(y, y0))
## All missing values
lx <- c(NA_real_, NA_real_)
y0 <- logSumExp_R(lx, na.rm = TRUE)
y <- logSumExp(lx, na.rm = TRUE)
print(y)
stopifnot(identical(y, -Inf))
stopifnot(all.equal(y, y0))
lx <- c(NA_real_, NA_real_)
y0 <- logSumExp_R(lx, na.rm = FALSE)
y <- logSumExp(lx, na.rm = FALSE)
print(y)
stopifnot(identical(y, NA_real_))
stopifnot(all.equal(y, y0))
## +Inf values
lx <- c(1, 2, +Inf)
y0 <- logSumExp_R(lx)
y <- logSumExp(lx)
print(y)
stopifnot(identical(y, +Inf))
stopifnot(all.equal(y, y0))
## First element is a missing value, cf. PR #33
lx <- c(NA_real_, 1)
y0 <- logSumExp_R(lx)
print(y0)
y <- logSumExp(lx, na.rm = FALSE)
print(y)
stopifnot(identical(y, NA_real_))
stopifnot(all.equal(y, y0))
y0 <- logSumExp_R(lx, na.rm = TRUE)
print(y0)
y <- logSumExp(lx, na.rm = TRUE)
print(y)
stopifnot(identical(y, 1))
stopifnot(all.equal(y, y0))
## Multiple -Inf values, cf. issue #84
lx <- c(-Inf, -Inf)
y0 <- logSumExp_R(lx)
y <- logSumExp(lx)
print(y)
stopifnot(identical(y, -Inf))
stopifnot(all.equal(y, y0))
lx <- c(-Inf, 5, -Inf)
y0 <- logSumExp_R(lx)
y <- logSumExp(lx)
print(y)
stopifnot(identical(y, 5))
stopifnot(all.equal(y, y0))
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.