segmentation fault in integrate (PR#812)

I tried to integrate numerically a function wich is similar to the
following:

dummy <- function(x) { exp(-1*x) * dnorm(x) }
dummy(-100)

[1] 0

dummy(-1000)

[1] NaN

dummy(-10000)

[1] NaN

If I choose the lower boundary to be too small integrate causes a
segmentation fault:

library(integrate)
integrate(dummy, -100, 0)$value

[1] 1.387143

integrate(dummy, -1000, 0)$value

[1] NaN
Warning message:
Ifail=1, maxpts was too small. Check the returned relerr! in: adapt(1,
lower, upper, minpts, maxpts, functn, eps, ...)

integrate(dummy, -10000, 0)$value

Segmentation fault

integrate(dummy, -1000000, 0)$value

dummy <- function(x) { a <-  exp(-1*x) * dnorm(x); ifelse(is.na(a), 0,

integrate(dummy, -100, 0)$value

integrate(dummy, -1000, 0, maxpts=1000)$value

I'm running R.1.2.0 on Debian Linux.

Regards
Achim

version

         _                
platform i686-pc-linux-gnu
arch     i686             
os       linux-gnu        
system   i686, linux-gnu  
status                    
major    1                
minor    2.0              
year     2000             
month    12               
day      15               
language R   



-------------------------

Achim Zeileis
Institut für Statistik
Technische Universität Wien

-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-devel mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !)  To: r-devel-request@stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._

I tried to integrate numerically a function wich is similar to the
following:

dummy <- function(x) { exp(-1*x) * dnorm(x) }
dummy(-100)

[1] 0

dummy(-1000)

[1] NaN

dummy(-10000)

[1] NaN

If I choose the lower boundary to be too small integrate causes a
segmentation fault:

library(integrate)
integrate(dummy, -100, 0)$value

[1] 1.387143

integrate(dummy, -1000, 0)$value

[1] NaN
Warning message:
Ifail=1, maxpts was too small. Check the returned relerr! in: adapt(1,
lower, upper, minpts, maxpts, functn, eps, ...)

integrate(dummy, -10000, 0)$value

Segmentation fault

Running Version 1.2.0 Patched (2001-01-07) and integrate Version  2.2-3
on SuSE 6.4 gives:

integrate(dummy, -1000000, 0)$value

[1] NaN
Warning message:
Ifail=1, maxpts was too small. Check the returned relerr! in: adapt(1,
lower, upper, minpts, maxpts, functn, eps, ...)

I has not able to reproduce the segfault ;-)

BTW: if one uses

dummy <- function(x) { a <-  exp(-1*x) * dnorm(x); ifelse(is.na(a), 0,

a) }

this returns

integrate(dummy, -100, 0)$value

[1] 1.387143

integrate(dummy, -1000, 0, maxpts=1000)$value

[1] 0.6935715

On Fri, 12 Jan 2001, Torsten Hothorn wrote:

I tried to integrate numerically a function wich is similar to the
following:

dummy <- function(x) { exp(-1*x) * dnorm(x) }
dummy(-100)

[1] 0

dummy(-1000)

[1] NaN

dummy(-10000)

[1] NaN

If I choose the lower boundary to be too small integrate causes a
segmentation fault:

library(integrate)
integrate(dummy, -100, 0)$value

[1] 1.387143

integrate(dummy, -1000, 0)$value

[1] NaN
Warning message:
Ifail=1, maxpts was too small. Check the returned relerr! in: adapt(1,
lower, upper, minpts, maxpts, functn, eps, ...)

integrate(dummy, -10000, 0)$value

Segmentation fault

Running Version 1.2.0 Patched (2001-01-07) and integrate Version  2.2-3
on SuSE 6.4 gives:

integrate(dummy, -1000000, 0)$value

[1] NaN
Warning message:
Ifail=1, maxpts was too small. Check the returned relerr! in: adapt(1,
lower, upper, minpts, maxpts, functn, eps, ...)

I has not able to reproduce the segfault ;-)

BTW: if one uses

dummy <- function(x) { a <-  exp(-1*x) * dnorm(x); ifelse(is.na(a), 0,

a) }

this returns

integrate(dummy, -100, 0)$value

[1] 1.387143

integrate(dummy, -1000, 0, maxpts=1000)$value

[1] 0.6935715

Better, try

dummy <- function(x) { exp(-1*x + dnorm(x, log=TRUE))}

that is, do your calculations on log scale as far as possible.  This is
more accurate than cop-ing out at Inf*0.

Even better, you can adjust the normal mean to cancel the linear term
exactly: dnorm(x, mu) has linear term -mu*x. So use

dummy <- function(x) exp(0.5)*dnorm(x, -1)

You will still need to adjust maxpts or find a sensible lower end.

(BTW, I don't get a segfault either, and we need you to debug it
on your system, I am afraid.)