R Programming (6) - Numerical integration

 

Numerical integration

Numerical integration

• It is frequently necessary to compute

  

• If we know the antiderivative F, then

  

• However for many function f, we don’t know the closed form of antiderivative.

• We approximate the integral on the divided subinterval:

  

Rectagular method

• Review the rectangular method

rect <- function(ftn, a, b, n = 100){
  h <- (b-a)/n
  x.vec <- seq(a, b, by = h)
  f.vec <- sapply(x.vec, ftn)
  h * sum(f.vec)
}

ftn6 <- function(x) return(4*x^3)
rect(ftn6, 0, 1, n=20)

## [1] 1.1025

Trapzoidal rule

• Approximating the are under \(y=f(x)\) over the subinterval \([x_i, x_{i+1}]\) by a trapezoid.

• The area of each trapezoid:

  

• Trapezoidal approximation:

  

trapezoid <- function(ftn, a, b, n = 100) {
  h <- (b-a)/n
  x.vec <- seq(a, b, by = h)
  f.vec <- sapply(x.vec, ftn)
  h*(f.vec[1]/2 + sum(f.vec[2:n]) + f.vec[n+1]/2)
}

ftn6 <- function(x) return(4*x^3)
trapezoid(ftn6, 0, 1, n=20)

## [1] 1.0025

Simpson’s rule

• Simpson’s rule subdivides the interval \([a,b]\) into \(n\) (even) subintervals

• and approximate \(f\) by a parabola (polynomial of degree 2).

• As an approximation to the area, we use

• Now assuming \(n\) is even, we add up the approximation for subintervals \([x_{2i}, x_{2i+2}]\) to obtain Simpson’s approximations \(S\).

• Notice that the \(f(x_i)\) for \(i\) odd are all weighted 4, while the \(f(x_i)\) for \(i\) even (except 0 and n) are weighted 2.

simpson_n <- function(ftn, a, b, n = 100) {
  n <- max(c(2*(n %/% 2), 4))
  h <- (b-a)/n
  x.vec1 <- seq(a+h, b-h, by = 2*h)
  x.vec2 <- seq(a+2*h, b-2*h, by = 2*h)
  f.vec1 <- sapply(x.vec1, ftn)
  f.vec2 <- sapply(x.vec2, ftn)
  h/3*(ftn(a) + ftn(b) + 4*sum(f.vec1) + 2*sum(f.vec2))   # return value
}

# or similarly
simpson_n <- function(ftn, a, b, n = 100) {
  n <- max(c(2*(n %/% 2), 4))
  h <- (b-a)/n
  
  x.vec <- seq(a, b, by = h)
  f.vec <- sapply(x.vec, ftn)
  
  h * 1/3 * (f.vec[1] + 4 * sum(f.vec[seq(2, n, by=2)]) + 2 * sum(f.vec[seq(3, n-1, by=2)]) + f.vec[n+1])
  
}

• Example

ftn6 <- function(x) return(4*x^3)
simpson_n(ftn6, 0, 1, 20)

## [1] 1

• Example

f <- function(x) {
  if (0 < x & x < 1) 1/simpson_n(function(x) exp(-x^3), 0, 1)*exp(-x^3)
  else 0
}

plot(seq(-1,2,0.01), sapply(seq(-1,2,0.01), f), xlab="x", ylab="f(x)", 'l')

• compute the mean

(m <- simpson_n(function(x) x*f(x), 0, 1))

## [1] 0.4317834

• compute the variance

simpson_n(function(x) (x-m)^2*f(x), 0, 1)

## [1] 0.0719255

• How to define the cumulative function F(x)?

pdf and cdf of the standard normal

Consider the probability density function and cumulative distribution of standard normal.

# probability density function
phi <- function(x) exp(-x^2 / 2) / sqrt(2*pi)

# cumulative distribution function
Phi <- function(z) {
    if (z < 0) 0.5 - simpson_n(phi, z, 0)
    else 0.5 + simpson_n(phi, 0, z)
}

z <- seq(-5, 5, by = 0.1)
phi.z <- sapply(z, phi)
Phi.z <- sapply(z, Phi)
plot(z, Phi.z, type  ="l", ylab="", main="phi(z) and Phi(z)")
lines(z, phi.z)

Functional programming

Return values are determined by f.vec:

• rectangular : h * sum(f.vec)
• trapzoid : h * (f.vec[1]/2 + sum(f.vec[2:n]) + f.vec[n+1]/2)
• Simpson : h * 1/3 * (f.vec[1] + 4 * sum(f.vec[seq(2, n, by=2)]) + 2 * sum(f.vec[seq(3, n-1, by=2)]) + f.vec[n+1])

Define each component as a function:

rect_method <- function(f.vec) sum(f.vec)

trapz_method <- function(f.vec) {
  n <- length(f.vec) - 1
  f.vec[1]/2 + sum(f.vec[2:n]) + f.vec[n+1]/2
}

simpson_method <- function(f.vec) {
  n <- length(f.vec) - 1
  (f.vec[1] +
    4 * sum(f.vec[seq(2, n, by=2)]) +
    2 * sum(f.vec[seq(3, n-1, by=2)]) +
    f.vec[n+1]) / 3
}

Now we can define general numerical integration function.

numerical_int <- function(ftn, a, b, n = 100, method){
  
  n <- max(c(2*(n %/% 2), 4))
  h <- (b-a)/n
  x.vec <- seq(a, b, by = h)
  f.vec <- sapply(x.vec, ftn)
  
  h * method(f.vec)
}

ftn6 <- function(x) return(4*x^3)

numerical_int(ftn6, 0, 1, method = rect_method)

## [1] 1.0201

numerical_int(ftn6, 0, 1, method = trapz_method)

## [1] 1.0001

numerical_int(ftn6, 0, 1, method = simpson_method)

## [1] 1

This programming tatic allow us to extend the numerical integrations as a new method is implemented.

For example, the Milne method:

milne_method <- function(f.vec) {
  n <- length(f.vec) - 1
  (4 * f.vec[1] +
    (-2) * sum(f.vec[seq(2, n, by=2)]) +
    8 * sum(f.vec[seq(3, n-1, by=2)]) +
    4 * f.vec[n+1]) * 1 / 3
}

numerical_int(ftn6, 0, 1, method = milne_method)

## [1] 1.0006

More about functional programming : http://adv-r.had.co.nz/Functional-programming.html