[R] integrate
Rafael A. Irizarry
ririzarr at jhsph.edu
Sun Jun 1 19:20:51 CEST 2003
Im tryng to understand an error i get with integrate. this is 1.7.0 on
solaris 2.8.
##i am trying to approximate an integral of this function,
f<-function(b) exp(-(b-mu)^2/(2*tau2))/(p-exp(b))*10^6
##with
tau2 <- .005;mu <- 7.96;p <- 2000
##from -inf to different upper limits. using
integrate(f,-Inf,log(p-exp(1)))
##i get the following error:
##Error in integrate(f, -Inf, log(p - exp(1))) :
## the integral is probably divergent
##whats confusing is that i only get this error when i use upper limit
##in the range [log(p-2),log(p-1)]
##try this:
x <- seq(1,10,.1)
y <- sapply(x,function(i){
x <- try(integrate(f,-Inf,log(p-i)),silent=TRUE)
if(class(x)=="try-error") return(NA) else return(x$val)
})
x[is.na(y)]
##if i use stop=FALSE, the curve is interpolated nicely
plot(x,y)
y <- sapply(x,function(i)
x <- integrate(f,-Inf,log(p-i),stop=FALSE)$val
)
lines(x,y)
##there doesnt appear to be anything strange about the function in
##the 2-3 region.
x <- seq(1,5,.1)
xx <- log(p-x)
plot(x,f(xx))
abline(v=c(2,3))
##any ideas on why this is happening?
##ps - same thing happens with:
##f <- function(b) exp(-(b-mu)^2/(2*tau2) - log(p-exp(b)))*10^6
##it doesnt with
##f<-function(b) exp(-(b-mu)^2/(2*tau2))/(p-exp(b))
More information about the R-help
mailing list