logspace.utils {statnet.common}

Utilities for performing calculations on logarithmic scale.

Description

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

Usage

log_sum_exp(logx, use_ldouble = FALSE)

log_mean_exp(logx, use_ldouble = FALSE)

lweighted.mean(x, logw)

lweighted.var(x, logw)

lweighted.cov(x, y, logw)

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.

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))

Author(s)

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))