Numerical Integration (Simpson’s Rule)

Abstract

Numerical integration is the study of how the numerical value of an integral can be found. Also called quadrature, which refers to nding a square whose area is the same as the area under a curve, it is one of the classical topics of numerical analysis. Of central interest is the process of approximating a deenite integral from values of the in-tegrand when exact mathematical integration is not available. The corresponding problem for multiple dimensional integration is known as multiple integration or cubature. Numerical integration has always been useful in bio-statistics to evaluate distribution functions and other quantities. Emphasis in recent years on Bayesian and empirical Bayesian methods and on mixture models has greatly increased the importance of numerical integration for computing likelihoods and posterior distributions and associated moments and derivatives. Many recent statistical methods are dependent especially on multiple integration , possibly in very high dimensions. This article describes classical quadrature methods and, more brieey, some of the more advanced methods for which software is widely available. The description of the elementary methods in this article borrows from introductory notes by Stewart 31]. An excellent general reference on numerical integration is 5]. More recent material can be found in 8] and 29]. Recent surveys of numerical integration with emphasis on statistical methods and applications are 10] and 9]. Trapezoidal Rule The simplest quadrature rule in wide use in the trapezoidal rule. Like many other methods, it has both a geometric and an analytic derivation. The idea of the geometric derivation is to approximate the area under the curve y = f(x) from x = a to x = b by the area of the trapezoid bounded by the points (a; 0), (b; 0), (a; f(a)) and (b; f(b)). This gives Z b a f(x)dx b ? a 2 ff(a) + f(b)g: The analytic derivation is to interpolate f(x) at a and b by a linear polynomial. The trapezoidal rule cannot be expected to give accurate results over a larger interval. However by summing the results of many applications of the trapezoidal rule over smaller intervals, we can obtain an accurate approximation to the integral over any interval. We begin by dividing a; b] into n equal intervals by the points a = x 0 < x 1 < < x n?1 < x n = b: Speciically, if h = b ? a n is the common length of …