You can divide your domain of integration into smaller intervals and then add up the individual contributions. This could improve the speed of adaptive Gauss-Kronrod quadrature used in integrate().
Ravi. ---------------------------------------------------------------------------- ------- Ravi Varadhan, Ph.D. Assistant Professor, The Center on Aging and Health Division of Geriatric Medicine and Gerontology Johns Hopkins University Ph: (410) 502-2619 Fax: (410) 614-9625 Email: [EMAIL PROTECTED] Webpage: http://www.jhsph.edu/agingandhealth/People/Faculty/Varadhan.html ---------------------------------------------------------------------------- -------- -----Original Message----- From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] On Behalf Of Santanu Pramanik Sent: Wednesday, August 22, 2007 6:41 PM To: [EMAIL PROTECTED] Cc: r-help@stat.math.ethz.ch; [EMAIL PROTECTED] Subject: Re: [R] integrate Hi Andy, Thank your very much for your input. I also tried something like that which gives a value close to 20, basically using the same trapezoidal rule. > sum(apply(as.matrix(seq(-10,10,by=0.1)),1,my.fcn))*0.1 [1] 20.17385 Actually my function is much more complicated than the posted example and it is taking a long time... Anyway, thanks again, Santanu JPSM, 1218J Lefrak Hall University of Maryland, College Park Phone: 301-314-9916 >>> <[EMAIL PROTECTED]> 8/22/2007 5:20:05 PM >>> As Duncan Murdoch mentioned in his reply, the problem is with the fact that your function is not really a properly defined function in the sense of "assigning a unique y to each x". The integrate function uses an adaptive quadrature routine which probably makes multiple calls to the function being integrated and expects to get the same y's for the same x's every time. If you want to get a number close to 20 (for your example) you need an integration routine which will use single evaluation of your "function" at each value of x. A simple method like rectangular approximation on a grid or the so-called trapezoidal rule will do just that. Here is a very crude prototype of such an integrator: integrate1 <- function(f, lower, upper){ f <- match.fun(f) xx <- seq(lower, upper, length=100) del <- xx[2] - xx[1] yy <- f(xx[-100]) return(del*sum(yy)) Now when you run integrate1(my.fun, -10, 10) you will get a number close to 20 but, of course, every time you do it you will get a different value. Hope this helps, Andy __________________________________ Andy Jaworski 518-1-01 Process Laboratory 3M Corporate Research Laboratory ----- E-mail: [EMAIL PROTECTED] Tel: (651) 733-6092 Fax: (651) 736-3122 "Santanu Pramanik" <[EMAIL PROTECTED] To .umd.edu> <r-help@stat.math.ethz.ch> Sent by: cc [EMAIL PROTECTED] at.math.ethz.ch Subject [R] integrate 08/22/2007 02:56 PM Hi, I am trying to integrate a function which is approximately constant over the range of the integration. The function is as follows: > my.fcn = function(mu){ + m = 1000 + z = 0 + z.mse = 0 + for(i in 1:m){ + z[i] = rnorm(1, mu, 1) + z.mse = z.mse + (z[i] - mu)^2 + } + return(z.mse/m) + } > my.fcn(-10) [1] 1.021711 > my.fcn(10) [1] 0.9995235 > my.fcn(-5) [1] 1.012727 > my.fcn(5) [1] 1.033595 > my.fcn(0) [1] 1.106282 > The function takes the value (approx) 1 over the range of mu. So the integration result should be close to 20 if we integrate over (-10, 10). But R gives: > integrate(my.fcn, -10, 10) 685.4941 with absolute error < 7.6e-12 > integrate(Vectorize(my.fcn), -10, 10) # this code never terminate I have seen in the "?integrate" page it is clearly written: If the function is approximately constant (in particular, zero) over nearly all its range it is possible that the result and error estimate may be seriously wrong. But this doesn't help in solving the problem. Thanks, Santanu JPSM, 1218J Lefrak Hall University of Maryland, College Park Phone: 301-314-9916 [[alternative HTML version deleted]] ______________________________________________ R-help@stat.math.ethz.ch 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. ______________________________________________ R-help@stat.math.ethz.ch 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. ______________________________________________ R-help@stat.math.ethz.ch 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.
