[R] Questions regarding "integrate" function

Le Wang ruser.le at gmail.com
Sat Nov 18 19:24:26 CET 2006


Hi there. Thanks for your time in advance.

I am using R 2.2.0 and OS: Windows XP.

My final goal is to calculate 1/2*integral of
(f1(x)^1/2-f2(x)^(1/2))^2dx (Latex codes:
$\frac{1}{2}\int^{{\infty}}_{\infty}
(\sqrt{f_1(x)}-\sqrt{f_2(x)})^2dx $.) where f1(x) and f2(x) are two
marginal densities.

My problem:

I have the following R codes using "adapt" package. Although "adapt"
function is mainly designed for more than 2 dimensions, the manual
says it will also call up "integrate" if the number of dimension
equals one. I feed in the data x1 and x2 and bandwidths h1 and h2.
These codes worked well when my final goal was to take double
integrals.

       integrand <- function(x) {
           # x input is evaluation point for x1 and x2, a 2x1 vector
           x1.eval <- x[1]
           x2.eval <- x[2]
           # n is the number of observations
           n <- length(x1)
           # x1 and x2 are the vectors read from data.dat
           # Compute the marginal densities
           f.x1 <- sum(dnorm((x1.eval-x1)/h1))/(n*h1)
           f.x2 <- sum(dnorm((x2.eval-x2)/h2))/(n*h2)
           # Return the integrand #
           return((sqrt(f.x1)-sqrt(f.x2))**2)

       }

       estimate<-0.5*adapt(1, lo=lo.default, up=up.default,
       minpts=minpts.default, maxpts=maxpts.default,
       functn=integrand, eps=eps.default, x1, x2,h1,h2)$value


But when I used it for one-dimension, it failed. Some of my
colleagues suggested getting rid of "x2.eval" in the "integrand"
because it is only one integral. But after I changed it, it still
didn't work. R gave the error msg: "evaluation of function gave a
result of wrong length"

I am not a frequent R user..although I looked up the mailing list
for a while and there were few postings asking similar questions, I
can't still figure out why my codes won't work. Any help will be
appreciated.

Le
-
~~~~~~~~~~~~~~~~~~~~~~
Le Wang, Ph.D
Population Center
University of Minnesota



More information about the R-help mailing list