[R] integrate

apjaworski at mmm.com apjaworski at mmm.com
Wed Aug 22 23:20:05 CEST 2007

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

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])

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 Jaworski
Process Laboratory
3M Corporate Research Laboratory
E-mail: apjaworski at mmm.com
Tel:  (651) 733-6092
Fax:  (651) 736-3122

             <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                                              

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.

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
PLEASE do read the posting guide
and provide commented, minimal, self-contained, reproducible code.

More information about the R-help mailing list