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 terminateI 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]]
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 terminateI 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]]
On 8/22/2007 3:54 PM, Santanu Pramanik wrote:> Hi, > I am trying to integrate a function which is approximately constant > over the range of the integration. The function is as follows:That's not a function of the input mu. It includes a random component: > my.fcn(10) [1] 0.9786558 > my.fcn(10) [1] 1.022467 You can't expect integrate() to return a sensible answer if you don't give it a function that returns consistent results. Duncan Murdoch> >> 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 at 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.
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>>> <apjaworski at mmm.com> 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: apjaworski at mmm.com Tel: (651) 733-6092 Fax: (651) 736-3122 "Santanu Pramanik" <spramanik at survey To .umd.edu> <r-help at stat.math.ethz.ch> Sent by: cc r-help-bounces at st 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 terminateI 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 at 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.