I want to measure the error in the estimation of a 2 dimensional density
function that is calculated using an integral but run into problems trying
to integrate with adaptIntegrate because the integrand also calls the
function adaptIntegrate. In particular I want
\int \hat{f}(x,y) - f(x,y) dx dy
where \hat{f}(x,y) = \int K(a,b, x, y) da db and in this simulation study I
know what the true value of f(x,y) is. I can use adaptIntegrate to
calculate \hat{f}(x1,y1) = adaptIntegrate( K, lowerLimit=c(a1,b1),
upperLimit=c(a2,b2), x=x1, y=y1)$integral which gives a real
number. However if I define the following integrand and try to use
adaptIntegrate on it I run into problems:
integrand <- function(x1,y1){
return(adaptIntegrate( K, lowerLimit=c(a1,b1), upperLimit=c(a2,b2),
x=x1, y=y1)$integral - f(x1,y1))
}
adaptIntegrate(integrand, lowerLimit = c(x1,y1), upperLimit = c(x2, y2))
Either adaptIntegrate evaluates the integrand at only one pair of (x,y)
values and returns the resulting integral or it decides it needs to look at
more than one point to achieve the desired accuracy but this seems to cause
an issue and it returns the following error:
REAL() can only be applied to a 'numeric', not a 'pairlist'
It seems like the adaptive part of the method is looking at what is inside
my integrand and running into problems when it tries to work out which
points to evaluate at. Although my integrand returns a real number when
called there might be an issue that adaptIntegrate first outputs a list and
then because I use adaptIntegrate(...)$integral I get the real number I
desire. Is there some way around that list format? Is that likely to cause
a problem in the integration?
Or any other deterministic multidimensional integrators? (I tried cuhre
from R2Cuba but it would give different results when I repeated the same
integration so I don´t want to use that).
Thanks for any suggestions!
[[alternative HTML version deleted]]
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
