Skip to content

Instantly share code, notes, and snippets.

@richarddmorey
Last active May 3, 2016 16:19
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
Star You must be signed in to star a gist
Save richarddmorey/3c77d0065983e31241bff3807482443e to your computer and use it in GitHub Desktop.
Representing real numbers using logarithms in R
###########
# S4 class to represent real numbers as logarithms
# Includes standard arithmetic operations +,-,*,/, and ^
# Richard D. Morey (richarddmorey@gmail.com)
# November 2014
###########
# Compute log(exp(a) + exp(b))
logExpAplusExpB <- function(a, b)
{
a + log1plusExpA( b - a )
}
# Compute log(1 + exp(a))
# From plogis.c
log1plusExpA <- function(a)
{
if( a <= 18 ) return( log1p( exp(a) ) )
if( a > 33.3 ) return( a )
# else: 18.0 < x <= 33.3 :
return( a + exp( -a ) )
}
# Compute log(exp(a) - exp(b))
logExpAminusExpB <- function(a, b)
{
a + pexp( a - b, log.p=TRUE )
}
setClass("logRepresentedReal",
representation(modulo = "numeric" ,
sign ="integer")
)
logRepresentedReal <- function( modulo, sign ){
new( "logRepresentedReal", modulo = modulo, sign = sign )
}
setValidity("logRepresentedReal", function(object){
if( abs(object@sign) != 1 ) return("Sign must be 1 or -1.")
return(TRUE)
})
setMethod("show", signature("logRepresentedReal"), function(object){
# Use the function in the BayesFactor package to print
# logarithmically-represented values nicely
require(BayesFactor, quietly = TRUE)
str = BayesFactor:::expString(object@modulo)
sgn = ""
if(object@sign < 0) sgn = "-"
cat(paste(" ", sgn, str, " ", sep=""))
})
setAs("logRepresentedReal", "numeric", function( from, to ){
as.double.logRepresentedReal(from)
})
as.double.logRepresentedReal <- function(x){
x@sign * exp(x@modulo)
}
setMethod("abs", signature("logRepresentedReal"), function(x){
x@sign = 1L
return(x)
})
setMethod("sign", signature("logRepresentedReal"), function(x){
x@sign
})
setMethod("-", signature("logRepresentedReal"), function(e1){
e1@sign = as.integer(-e1@sign)
return(e1)
})
setMethod('>', signature("logRepresentedReal", "logRepresentedReal"), function(e1, e2){
if( (e1@modulo == e2@modulo) & (e1@sign == e2@sign) ) return(FALSE)
if(sign(e1) > sign(e2)) return(TRUE)
if(sign(e1) < sign(e2)) return(FALSE)
if(sign(e1) == 1) return(e1@modulo > e2@modulo)
if(sign(e1) == -1) return(e1@modulo < e2@modulo)
})
setMethod('<', signature("logRepresentedReal", "logRepresentedReal"), function(e1, e2){
if( (e1@modulo == e2@modulo) & (e1@sign == e2@sign) ) return(FALSE)
return(!(e1>e2))
})
setGeneric("is.zero", function(x){
standardGeneric("is.zero")
})
setMethod("is.zero", signature("logRepresentedReal"), function(x){
if( !is.infinite(x@modulo) ) return(FALSE)
if(sign(x) != sign(x@modulo)) return(TRUE)
return(FALSE)
})
setMethod("is.infinite", signature("logRepresentedReal"), function(x){
if( !is.infinite(x@modulo) ) return(FALSE)
if(sign(x) == sign(x@modulo)) return(TRUE)
return(FALSE)
})
setMethod("is.finite", signature("logRepresentedReal"), function(x){
return(!is.infinite(x))
})
setMethod('+', signature("logRepresentedReal", "logRepresentedReal"), function(e1, e2){
if(e1@sign == 1 & e2@sign == -1 ) return(e1 - -e2)
if(e1@sign == -1 & e2@sign == -1 ) return( -(-e1 + -e2) )
if(e1@sign == -1 & e2@sign == 1 ) return(e2 - -e1)
logsum = logExpAplusExpB(e1@modulo, e2@modulo)
logRepresentedReal( logsum, 1L )
})
setMethod('+', signature("logRepresentedReal", "numeric"), function(e1, e2){
if( e2 == 0 ) return(e1)
e2 = logRepresentedReal( log(abs(e2)), as.integer(sign(e2)) )
e1 + e2
})
setMethod('+', signature("numeric","logRepresentedReal"), function(e1, e2){
e2 + e1
})
setMethod('-', signature("logRepresentedReal", "logRepresentedReal"), function(e1, e2){
if(e1@sign == 1 & e2@sign == -1 ) return( e1 + -e2 )
if(e1@sign == -1 & e2@sign == -1 ) return( -abs(e1) + abs(e2) )
if(e1@sign == -1 & e2@sign == 1 ) return( -(e2 + -e1) )
if(e1 > e2){
logdiff = logExpAminusExpB(e1@modulo, e2@modulo)
return(logRepresentedReal( logdiff, 1L ))
}else if(e1 < e2){
logdiff = logExpAminusExpB(e2@modulo, e1@modulo)
return(logRepresentedReal( logdiff, -1L ))
}else{
return(logRepresentedReal( -Inf, 1L ))
}
})
setMethod('-', signature("logRepresentedReal", "numeric"), function(e1, e2){
if( e2 == 0 ) return(e1)
e2 = logRepresentedReal( log(abs(e2)), as.integer(sign(e2)) )
e1 - e2
})
setMethod('-', signature("numeric","logRepresentedReal"), function(e1, e2){
-e2 + e1
})
setMethod('*', signature("numeric", "logRepresentedReal"), function(e1, e2){
sgn = sign(e1) * sign(e2)
modu = log(abs(e1)) + e2@modulo
logRepresentedReal( modu, as.integer(sgn) )
})
setMethod('*', signature("logRepresentedReal", "numeric"), function(e1, e2){
e2 * e1
})
setMethod('*', signature("logRepresentedReal", "logRepresentedReal"), function(e1, e2){
sgn = sign(e1) * sign(e2)
modu = e1@modulo + e2@modulo
logRepresentedReal( modu, as.integer(sgn) )
})
setMethod('/', signature("logRepresentedReal", "logRepresentedReal"), function(e1, e2){
sgn = sign(e1) * sign(e2)
modu = e1@modulo - e2@modulo
logRepresentedReal( modu, as.integer(sgn) )
})
setMethod('/', signature("numeric", "logRepresentedReal"), function(e1, e2){
if(e1 == 1){
e2@modulo = -e2@modulo
return(e2)
}else{
return(e1 * (1/e2))
}
})
setMethod('/', signature("logRepresentedReal", "numeric"), function(e1, e2){
1/(e2/e1)
})
setMethod('^', signature("logRepresentedReal", "numeric"), function(e1, e2){
if( e2 == 0 )
return( logRepresentedReal( 0, 1L ) )
if( sign(e1)==1 | e2%%2 == 0 ){
sgn = 1L
}else{
sgn = -1L
}
if(e2 < 0){
e2 = abs(e2)
e1 = 1/e1
}
modu = e1@modulo * e2
logRepresentedReal( modu, as.integer(sgn) )
})
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment