home *** CD-ROM | disk | FTP | other *** search
/ Big Green CD 8 / BGCD_8_Dev.iso / NEXTSTEP / UNIX / Educational / R-0.49-MI / R-0.49-I / demos / language / recursion next >
Encoding:
Text File  |  1996-11-24  |  2.0 KB  |  81 lines
  1. # Adaptive integration:  Venables and Ripley pp. 105-110
  2. # This is the basic integrator.
  3.  
  4. area <- function(f, a, b, ..., fa = f(a, ...), fb = f(b, ...), limit
  5.      = 10, eps = 1.e-5)
  6. {
  7.     h <- b - a      
  8.     d <- (a + b)/2
  9.     fd <- f(d, ...) 
  10.     a1 <- ((fa + fb) * h)/2
  11.     a2 <- ((fa + 4 * fd + fb) * h)/6
  12.     if(abs(a1 - a2) < eps)  
  13.         return(a2)  
  14.     if(limit == 0) {
  15.         warning(paste("iteration limit reached near x = ",
  16.             d))
  17.         return(a2)
  18.     }       
  19.     area(f, a, d, ..., fa = fa, fb = fd, limit = limit - 1,
  20.         eps = eps) + area(f, d, b, ..., fa = fd, fb =
  21.         fb, limit = limit - 1, eps = eps)
  22. }
  23.  
  24.  
  25. # The function to be integrated
  26.  
  27. fbeta <- function(x, alpha, beta)
  28. {
  29.     x^(alpha - 1) * (1 - x)^(beta - 1)
  30. }
  31.  
  32.  
  33. # Compute the approximate integral, the exact integral and the error
  34.  
  35. b0 <- area(fbeta, 0, 1, alpha=3.5, beta=1.5)
  36. b1 <- exp(lgamma(3.5) + lgamma(1.5) - lgamma(5))
  37. c(b0, b1, b0-b1)
  38.  
  39.  
  40. # Modify the function so that it records where it was evaluated
  41.  
  42. fbeta.tmp <- function (x, alpha, beta) 
  43. {
  44.     val <<- c(val, x)
  45.     x^(alpha - 1) * (1 - x)^(beta - 1)
  46. }
  47.  
  48.  
  49. # Recompute and plot the evaluation points.
  50.  
  51. val <- NULL
  52. b0 <- area(fbeta.tmp, 0, 1, alpha=3.5, beta=1.5)
  53. plot(val, fbeta(val, 3.5, 1.5), pch=0)
  54.  
  55.  
  56. # Better programming style -- renaming the function will have no effect.
  57. # The use of "Recall" as in V+R is VERY black magic.  You can get the
  58. # same effect transparently by supplying a wrapper function.
  59. # This is the approved Abelson+Sussman method.
  60.  
  61. area <- function(f, a, b, ..., limit=10, eps=1e-5) {
  62.     area2 <- function(f, a, b, ..., fa = f(a, ...), fb = f(b, ...),
  63.             limit = limit, eps = eps) {
  64.         h <- b - a
  65.         d <- (a + b)/2
  66.         fd <- f(d, ...) 
  67.         a1 <- ((fa + fb) * h)/2
  68.         a2 <- ((fa + 4 * fd + fb) * h)/6
  69.         if(abs(a1 - a2) < eps)
  70.             return(a2)
  71.         if(limit == 0) {
  72.             warning(paste("iteration limit reached near x =", d))
  73.             return(a2)
  74.         }
  75.         area2(f, a, d, ..., fa = fa, fb = fd, limit = limit - 1,
  76.             eps = eps) + area2(f, d, b, ..., fa = fd, fb =
  77.             fb, limit = limit - 1, eps = eps)
  78.     }
  79.     area2(f, a, b, ..., limit=limit, eps=eps)
  80. }
  81.