Created
April 25, 2011 20:15
-
-
Save doobwa/941125 to your computer and use it in GitHub Desktop.
Computing log-sum-exponential
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
log.sum.exp<- function(x) { | |
# Computes log(sum(exp(x)) | |
# Uses offset trick to avoid numeric overflow: http://jblevins.org/notes/log-sum-exp | |
offset <- x[which(x == max(abs(x)))] | |
log(sum(exp(x - offset))) + offset | |
} |
As mentioned in http://jblevins.org/notes/log-sum-exp and http://scicomp.stackexchange.com/questions/1122/how-to-add-large-exponential-terms-reliably-without-overflow-errors, one can use the following code:
log.sum.exp<- function(x) {
# Computes log(sum(exp(x))
# Uses offset trick to avoid numeric overflow: http://jblevins.org/notes/log-sum-exp
if ( max(abs(x)) > max(x) )
offset <- min(x)
else
offset <- max(x)
log(sum(exp(x - offset))) + offset
}
And with the possibility to apply only over certain margins/axis:
log.sum.exp <- function(x, margin=NULL) {
# Computes log(sum(exp(x))
# Uses offset trick to avoid numeric overflow: http://jblevins.org/notes/log-sum-exp
innerf <- function(x) {
if ( max(abs(x)) > max(x) )
offset <- min(x)
else
offset <- max(x)
log(sum(exp(x - offset))) + offset
}
if (is.null(margin)) {
return(innerf(x))
}
else {
return(apply(x, margin, innerf))
}
}
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Hi,
Thank you for the snippet, but this snippet doesnt work for numerical underflow problems.
Regards
Måns