Rectangular Integration (a.k.a. The Midpoint Rule)

Posted on January 20, 2014 by Eric Cai - The Chemical Statistician in R bloggers | 0 Comments

[This article was first published on The Chemical Statistician » R programming, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)

Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

Introduction

Continuing on the recently born series on numerical integration, this post will introduce rectangular integration. I will describe the concept behind rectangular integration, show a function in R for how to do it, and use it to check that the $Beta(2, 5)$ distribution actually integrates to 1 over its support set. This post follows from my previous post on trapezoidal integration.

Image courtesy of Qef from Wikimedia Commons.

Conceptual Background of Rectangular Integration (a.k.a. The Midpoint Rule)

Rectangular integration is a numerical integration technique that approximates the integral of a function with a rectangle. It uses rectangles to approximate the area under the curve. Here are its features:

The rectangle’s width is determined by the interval of integration.
- One rectangle could span the width of the interval of integration and approximate the entire integral.
- Alternatively, the interval of integration could be sub-divided into $n$ smaller intervals of equal lengths, and $n$ rectangles would used to approximate the integral; each smaller rectangle has the width of the smaller interval.
The rectangle’s height is the function’s value at the midpoint of its base.
Within a fixed interval of integration, the approximation becomes more accurate as more rectangles are used; each rectangle becomes narrower, and the height of the rectangle better captures the values of the function within that interval.

This method is inspired by the use of Riemann sums to calculate the integral. Roughly speaking, the limit of the Riemann sums of a function as partitions become finer is the Riemann integral. A function is Riemann-integrable if this limit exists, and the Riemann sum becomes closer to the Riemann integral with a sufficiently fine partition.

Recall that the area of a rectangle is just the base multiplied by the height.

$\text{Area of Rectangle} = \text{base} \times \text{height}$ .

To adapt this to a rectangle that approximates a function $f(x)$ over an interval $[a, b]$ ,

$\text{Area of Rectangle} = (b - a) \times f[0.5(b - a)]$

To approximate the integral of $f(x)$ with $n$ rectangles over the interval $[a, b]$ ,

$\int_{a}^{b} f(x)\, dx \approx \sum_{i=1}^{n-1}(x_{i+1} - x_{i}) \times f[0.5(x_i + x_{i+1})]$

R Script for Integrating $Beta(2, 5)$ Over $[0, 1]$ .

The following code implements rectangular integration to check that the $Beta(2, 5)$ PDF actually integrates to 1 over its support set. It defines 2 functions, creates a vector of the support set, and finally provides the integral. The 2 functions are:

dbeta.2.5(): This produces the PDF values for the Beta(2, 5) distribution. It is needed to calculate the midpoints.
rectangular.integration(): This implements rectangular integration and numerically approximates the integral. dbeta.2.5() is eventually fed into rectangular.integration as its second argument.

After the 2 functions are defined, the support set is created in beta.support, and the integral is computed.

Notice that the function itself is needed as an input for rectangular.integration(); the values of the function on the support set are not given as an input (in contrast to my example for trapezoidal integration. Thus, rectangular.integration() can only be implemented if the function values at the midpoints can be obtained. This may not always be possible, especially when some new function without an analytical expression needs to be integrated and only its values at certain points along its input are known. In my example, I am integrating the PDF of the Beta(2, 5) distribution, and any value of the PDF can be calculated in R, so this works for rectangular.integration().

If your integrand cannot be evaluated at the midpoints of your intervals, you can modify rectangular integration to use the function’s value at either the left boundary or right boundary as the rectangle’s height. These 2 values can be taken directly from the function values, but the resulting approximation is not as good as using the midpoint rule. Nonetheless, rectangular.integration() can be modified to do so, and I will leave this as an exercise for you to do.

If you try to “guess” $f[(a+b)/2]$ by taking the average of $f(a)$ and $f(b)$ (i.e. using $0.5 \times [f(a) + f(b)]$ ), then the resulting method is just trapezoidal integration.

##### Rectangular Integration (a.k.a. Midpoint Rule)
##### By Eric Cai - The Chemical Statistician
# clear all variables
rm(list = ls(all.names = TRUE))

### define a function that produces the PDF values for the Beta(2, 5) distribution
### Input: the quantiles
dbeta.2.5 = function(quantiles)
{
 # check if input is numeric
 if (!is.numeric(quantiles))
 {
 stop('The argument is not numeric.')
 }

 return(dbeta(quantiles, 2, 5))
}

### Define a function for implementing rectangular integration
### Input #1: x - the points along the interval of integration
### Input #2: f - the function
### Output: the integral as approximated by rectangular integration (a.k.a. midpoint rule)
rectangular.integration = function(x, f)
{
 # check if the variable of integration is numeric
 if (!is.numeric(x))
 {
 stop('The first argument is not numeric.')
 }

 # check if f is a function
 if (!is.function(f))
 {
 stop('The second argument is not a function.')
 }

 ### finish checks

 # obtain length of variable of integration and integrand
 n.points = length(x)

 # midpoints
 midpoints = 0.5*(x[2:n.points] + x[1:(n.points-1)])

 # function evaluated at midpoints
 f.midpoints = f(midpoints)

 # calculate the widths of the intervals between adjacent pairs of points along the variable of integration
 interval.widths = x[2:n.points] - x[1:(n.points-1)]

 # implement rectangular integration
 # calculate the sum of all areas of rectangles that are used to approximate the integral
 rectangular.integral = sum(interval.widths * f.midpoints)

 # print the definite integral
 return(rectangular.integral)
}

# points along support set for Beta(2, 5)
beta.support = seq(0, 1, by = 0.005)

# calculate integral of Beta(2, 5) PDF over its support set
rectangular.integration(beta.support, dbeta.2.5)
[1] 1.000031

Just like trapezoidal integration, the rectangular integral is not exactly 1. This error is caused by

the finite number of points that I used to define the Beta(2, 5) PDF, which leads to the finite number of rectangles that can be used to approximate the integral
the inherent error in rectangular integration

As mentioned previously, I will discuss the error analysis of trapezoidal integration and other numerical integration techniques in a later post.