Utilities for performing calculations on logarithmic scale.

Source: R/logspace.utils.R

A small suite of functions to compute sums, means, and weighted means on logarithmic scale, minimizing loss of precision.

log_sum_exp(logx, use_ldouble = FALSE)

log_mean_exp(logx, use_ldouble = FALSE)

lweighted.mean(x, logw)

lweighted.var(x, logw, onerow = NA)

lweighted.cov(x, y, logw, onerow = NA)

log1mexp(x)

Arguments

logx

Numeric vector of \(\log(x)\), the natural logarithms of the values to be summed or averaged.

use_ldouble

Whether to use long double precision in the calculation. If TRUE, 's C built-in logspace_sum() is used. If FALSE, the package's own implementation based on it is used, using double precision, which is (on most systems) several times faster, at the cost of precision.

x, y

Numeric vectors or matrices of \(x\) and \(y\), the (raw) values to be summed, averaged, or whose variances and covariances are to be calculated.

logw

Numeric vector of \(\log(w)\), the natural logarithms of the weights.

onerow

If given a matrix or matrices with only one row (i.e., sample size 1), var() and cov() will return NA. But, since weighted matrices are often a product of compression, the same could be interpreted as a variance of variables that do not vary, i.e., 0. This argument controls what value should be returned.

Value

The functions return the equivalents of the R expressions given below, but faster and with less loss of precision.

Functions

  • log_sum_exp(): log(sum(exp(logx)))

  • log_mean_exp(): log(mean(exp(logx)))

  • lweighted.mean(): weighted mean of x: sum(x*exp(logw))/sum(exp(logw)) for x scalar and colSums(x*exp(logw))/sum(exp(logw)) for x matrix

  • lweighted.var(): weighted variance of x: crossprod(x-lweighted.mean(x,logw)*exp(logw/2))/sum(exp(logw))

  • lweighted.cov(): weighted covariance between x and y: crossprod(x-lweighted.mean(x,logw)*exp(logw/2), y-lweighted.mean(y,logw)*exp(logw/2))/sum(exp(logw))

  • log1mexp(): log(1-exp(-x)) for x >= 0 (a wrapper for the eponymous C macro provided by R)

Author

Pavel N. Krivitsky

Examples

x <- rnorm(1000)
stopifnot(all.equal(log_sum_exp(x), log(sum(exp(x))), check.attributes=FALSE))
stopifnot(all.equal(log_mean_exp(x), log(mean(exp(x))), check.attributes=FALSE))

logw <- rnorm(1000)
stopifnot(all.equal(m <- sum(x*exp(logw))/sum(exp(logw)),lweighted.mean(x, logw)))
stopifnot(all.equal(sum((x-m)^2*exp(logw))/sum(exp(logw)),
                    lweighted.var(x, logw), check.attributes=FALSE))

x <- cbind(x, rnorm(1000))
stopifnot(all.equal(mx <- colSums(x*exp(logw))/sum(exp(logw)),
                    lweighted.mean(x, logw), check.attributes=FALSE))
stopifnot(all.equal(crossprod(t(t(x)-mx)*exp(logw/2))/sum(exp(logw)),
                    lweighted.var(x, logw), check.attributes=FALSE))


y <- cbind(x, rnorm(1000))
my <- colSums(y*exp(logw))/sum(exp(logw))
stopifnot(all.equal(crossprod(t(t(x)-mx)*exp(logw/2), t(t(y)-my)*exp(logw/2))/sum(exp(logw)),
                    lweighted.cov(x, y, logw), check.attributes=FALSE))
stopifnot(all.equal(crossprod(t(t(y)-my)*exp(logw/2), t(t(x)-mx)*exp(logw/2))/sum(exp(logw)),
                    lweighted.cov(y, x, logw), check.attributes=FALSE))

x <- rexp(1000)
stopifnot(isTRUE(all.equal(log1mexp(x), log(1-exp(-x)))))