integrate function fails! (PR#1718)

[Thomas Lumley]

On Fri, 28 Jun 2002 jechacon@unex.es wrote:

> Full_Name: José Enrique Chacón
> Version: 1.5.0 and 1.3.1
> OS: Windows Millenium
> Submission from: (NULL) (158.49.28.155)
>
>
> Dear reader:
> I was trying to evaluate the L2 error produced when estimating the density
> function N(0,1) from a sample of size 100 using a kernel density estimate. It
> produced a strange value. You can reproduce the process by typing
>
> samp<-rnorm(100)
> fker<-function(x,h){sum(dnorm((rep(x,100)-samp)/h))/(100*h)}
> integrand<-function(x){(fker(x,5)-dnorm(x))^2}
> err<-integrate(integrand,-Inf,Inf)
>
> This value err of the L2 error is generally far from accurate. fker is the
> kernel density estimate with bandwidth h and Gaussian kernel. Obviously, h=5 is
> a very large choice far from the optimum which is around 0.4, but it makes it
> more clear to see the failure. You just have to evaluate integrand(0); then, you
> know that the integral between -0.5 and 0.5 cannot be larger than integrand(0),
> which was usually around 0.1. But the value that integrate gives is around 1.6!
> Anyway, there is a library called rmutil (I don't remember where I got it, but I
> found it searching in the forums) which contains a function called int that
> works well.

This may illustrate what's going on...

> x<-seq(-1,1,length=10)
> sapply(x,integrand)
 [1] 0.02741093 0.04742150 0.06987683 0.08968701 0.10136155 0.10132284
 [7] 0.08957798 0.06971705 0.04723830 0.02723325
> integrand(x)
 [1] 0.2865455 0.2327631 0.1895539 0.1599095 0.1449956 0.1449956 0.1599095
 [8] 0.1895539 0.2327631 0.2865455


If you look a the help page for integrate you see

       f: An R function taking a numeric first argument and returning a
          numeric vector the same length. Returning a non-finite
          element will generate an error.

so the function must be properly vectorised for integrate to work. Your
function isn't.  There is a check in integrate() that the return value is
the right length, but there's really no way to check that the value is
correct.

If you make the function vectorised it does seem to work
> correct.integrand<-function(x) sapply(x,integrand)
> integrate(correct.integrand,-Inf,Inf)
0.1839647 with absolute error < 8.9e-08


	-thomas

Thomas Lumley			Asst. Professor, Biostatistics
tlumley@u.washington.edu	University of Washington, Seattle


-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._