[R] integrate

[apjaworski at mmm.com]

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 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.