logspace.utils | R Documentation |
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)
logx |
Numeric vector of |
use_ldouble |
Whether to use |
x , y |
Numeric vectors or matrices of |
logw |
Numeric vector of |
onerow |
If given a matrix or matrices with only one row
(i.e., sample size 1), |
The functions return the equivalents of the R expressions given below, but faster and with less loss of precision.
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))
Pavel N. Krivitsky
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))
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.