[R] integrate
Ravi Varadhan
rvaradhan at jhmi.edu
Thu Aug 23 00:48:38 CEST 2007
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: rvaradhan at jhmi.edu
Webpage: http://www.jhsph.edu/agingandhealth/People/Faculty/Varadhan.html
----------------------------------------------------------------------------
--------
-----Original Message-----
From: r-help-bounces at stat.math.ethz.ch
[mailto:r-help-bounces at stat.math.ethz.ch] On Behalf Of Santanu Pramanik
Sent: Wednesday, August 22, 2007 6:41 PM
To: apjaworski at mmm.com
Cc: r-help at stat.math.ethz.ch; r-help-bounces at stat.math.ethz.ch
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
>>> <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.
______________________________________________
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.
More information about the R-help
mailing list