--------------------------------------------------------------
Laurent Gautier CBS, Building 208, DTU
PhD. Student DK-2800 Lyngby,Denmark
tel: +45 45 25 24 89 http://www.cbs.dtu.dk/laurent
On Wed, Aug 07, 2002 at 12:28:28AM +0900, Lu Chi-Hsien Joseph wrote:
> Dear R-users,
>
> The numerical integration example given in Gentleman and Ihaka (2000),
> "Lexical Scope and Statistical Computing," JCGS, 9, 491-508,
> is very interesting and helpful in understanding how lexical scope
> is about.
> However, I got some questions that I just can't figure out.
> First all, allow me to copy the two functions given by the authors:
>
> midpoint <- function(f, a, b, n=100) {
> h <- (b - a)/n
> (b - a)*mean(sapply(seq(a + h/2, b - h/2, length=n), f))
> }
>
> bv.integrate <- function(f, a, b, n=100, rule=midpoint) {
> g <- function(y) {
> fx <- function(x) f(x, y)
> rule(fx, a[2], b[2], n)
> }
> rule(g, a[1], b[1])
> }
>
> I modified the function name from "integrate" to "bv.integrate".
> To integrate f(x,y)=sin(x + y) over the unit square,
> it works:
> > bv.integrate(function(x, y) sin(x + y), c(0,0), c(1,1))
> [1] 0.773651
> (compare with .7736445 calculated from 2sin(1)*(1 - cos(1)))
>
> Then I write a function using trapezoidal rule as follows:
>
> trapezoidal <- function(f, a, b, n=100) {
> h <- (b - a)/n
> i <- 1:n
> sum(f(a + (i - 1)*h), f(a + i*h))*h/2
> }
>
> I checked with
> > trapezoidal(sin, 0, 1)
> [1] 0.4596939
> > 1 - cos(1)
> [1] 0.4596977
>
> and
> > trapezoidal(dnorm, -3, 3)
> [1] 0.9972922
> > diff(pnorm(c(-3, 3)))
> [1] 0.9973002
>
> Then, this is what I got by plugged in "trapezoidal" for the "rule"
> in bv.integrate():
> > bv.integrate(function(x, y) sin(x + y), c(0,0), c(1,1), rule=trapezoidal)
> [1] 0.007080675
>
> Hundred time smaller!
> Answer can be "improved" by increasing n, but always 100 times smaller.
>
> Anything wrong with the way I wrote trapezoidal()?
> I just can't figure out any difference, in terms of input/output,
> between midpoint() and trapezoidal().
>
> Then I tried to use integrate() as the rule in bv.integrate(),
> by defining a wrap-up function:
>
> intg <- function(f, a, b, n=100) integrate(f, a, b, n)$value
>
> to have the value as the only output.
> But, I got
> > bv.integrate(function(x, y) sin(x + y), c(0,0), c(1,1), rule=intg)
> Error in integrate(f, a, b, n) : evaluation of function gave a result
> of wrong length
>
> I tried to use debugger() to find out what's wrong,
> but not succeeded.
> There must be something I missed in understanding the lexical scope,
> however, what's going on with the function I wrote?
> What should be the way to correct them?
>
> Certainly, I know integrate() should be the one to use.
> These are questions from my preparation of programming exercise
> of using R in my statistical computation class.
>
> I'll be very appreciated for anyone who might help me on these
> questions.
>
> Best regards,
>
> C. Joseph Lu
> Department of Statistics
> National Cheng-Kung University
> Tainan, Taiwan, ROC
> -.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
> r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
> Send "info", "help", or "[un]subscribe"
> (in the "body", not the subject !) To: r-help-request at stat.math.ethz.ch
> _._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !) To: r-help-request at stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._