[R] integrate

[Santanu Pramanik]

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