7. Quadrature and Interpolation #

Contents

Quadrature and Interpolation

7.1. Overview #

Julia has both a large number of useful, well written libraries and many incomplete poorly maintained proofs of concept.

A major advantage of Julia libraries is that, because Julia itself is sufficiently fast, there is less need to mix in low level languages like C and Fortran.

As a result, most Julia libraries are written exclusively in Julia.

Not only does this make the libraries more portable, it makes them much easier to dive into, read, learn from and modify.

See general, data, and statistical packages and optimization, solver, and related packages for more domain specific packages.

In this section we will explore the related concepts of Quadrature and Interpolation.

using LinearAlgebra, Statistics, Distributions
using QuadGK, FastGaussQuadrature, SpecialFunctions
using Interpolations, Plots

7.2. Numerical Integration #

Many applications require directly calculating a numerical derivative and calculating expectations.

7.2.1. Adaptive Quadrature#

A high accuracy solution for calculating numerical integrals is QuadGK.

using QuadGK
@show value, tol = quadgk(cos, -2π, 2π);

(value, tol) = quadgk(cos, -2π, 2π) = (-1.8403051655192592e-17, 1.0683776742236108e-23)

This is an adaptive Gauss-Kronrod integration technique that’s relatively accurate for smooth functions.

However, its adaptive implementation makes it slow and not well suited to inner loops.

7.2.2. Gauss Legendre#

Alternatively, many integrals can be done efficiently with (non-adaptive) Gaussian quadrature.

For example, using FastGaussQuadrature.jl

x, w = FastGaussQuadrature.gausslegendre(100_000); # i.e. find 100,000 nodes

# integrates f(x) = x^2 from -1 to 1
f(x) = x^2
@show dot(w, f.(x)); # calculate integral

dot(w, f.(x)) = 0.6666666666666665

With the FastGaussQuadrature package you may need to deal with affine transformations to the non-default domains yourself.

7.2.3. Gauss Legendre#

Another commonly used quadrature well suited to random variables with bounded support is Gauss–Jacobi quadrature.

It provides nodes $s_{n} \in [- 1, 1]$ and weights $ω_{n}$ for

\int_{- 1}^{1} g (s) (1 - s)^{a} (1 + s)^{b} d s \approx \sum_{n = 1}^{N} ω_{n} g (s_{n}) .

For $X \sim Beta (α, β)$ ,

E [f (X)] = \int_{0}^{1} f (x) \frac{x^{α - 1} (1 - x)^{β - 1}}{B (α, β)} d x,

with the change of variables $s = 2 x - 1$ (so $x = (s + 1) / 2$ , $d x = d s / 2$ ). This yields Gauss–Jacobi exponents $a = β - 1$ , $b = α - 1$ and a factor $C = 2^{- (α + β - 1)} / B (α, β)$ :

E [f (X)] \approx C \sum_{n = 1}^{N} ω_{n} f (\frac{s_{n} + 1}{2}) .

function gauss_jacobi(F::Beta, N)
    s, wj = FastGaussQuadrature.gaussjacobi(N, F.β - 1, F.α - 1)
    x = (s .+ 1) ./ 2
    C = 2.0^(-(F.α + F.β - 1.0)) / SpecialFunctions.beta(F.α, F.β)
    w = C .* wj
    return x, w
end

N = 32
F = Beta(2, 2)
F2 = Beta(0.5, 1.2)
x, w = gauss_jacobi(F, N)
x2, w2 = gauss_jacobi(F2, N)
f(x) = x^2
@show dot(w, f.(x)), mean(f.(rand(F, 1000)))
@show dot(w2, f.(x2)), mean(f.(rand(F2, 1000)));

(dot(w, f.(x)), mean(f.(rand(F, 1000)))) = (0.29999999999999993, 0.29447371629138674)

(dot(w2, f.(x2)), mean(f.(rand(F2, 1000)))) = (0.16339869281045832, 0.177096704880967)

Many expectations are of the form $E [f (X)] \approx \sum_{n = 1}^{N} w_{n} f (x_{n})$ where $X \sim N (0, 1)$ . Alternatively, integrals of the form $\int_{- \infty}^{\infty} f (x) e x p (- x^{2}) d x$ .

Gauss-hermite quadrature provides weights and nodes of this form, where the normalize = true argument provides the appropriate rescaling for the normal distribution.

Through a change of variables this can be used to calculate expectations of $N (μ, σ^{2})$ variables, through

function gauss_hermite_normal(N::Integer, mu::Real, sigma::Real)
  s, w = FastGaussQuadrature.gausshermite(N; normalize = true)

  # X ~ N(mu, sigma^2), X \sim mu + sigma N(0,1) we transform the standard-normal nodes
  x = mu .+ sigma .* s
  return x, w
end

N = 32
x, w = gauss_hermite_normal(N, 1.0, 0.1)
x2, w2 = gauss_hermite_normal(N, 0.0, 0.05)
f(x) = x^2
@show dot(w, f.(x)), mean(f.(rand(Normal(1.0, 0.1), 1000)))
@show dot(w2, f.(x2)), mean(f.(rand(Normal(0.0, 0.05), 1000))); 

(dot(w, f.(x)), mean(f.(rand(Normal(1.0, 0.1), 1000)))) = (1.01, 1.0119555974422323)

(dot(w2, f.(x2)), mean(f.(rand(Normal(0.0, 0.05), 1000)))) = (0.0024999999999999996, 0.0025550146158386127)

7.4. Gauss Laguerre #

Another quadrature scheme appropriate integrals defined on $[0, \infty]$ is Gauss Laguerre, which approximates integrals of the form $\int_{0}^{\infty} f (x) x^{α} \exp (- x) d x \approx \sum_{n = 1}^{N} w_{n} f (x_{n})$ .

One application is to calculate expectations of exponential variables. The PDF of an exponential distribution with parameter $θ$ is $f (x; θ) = \frac{1}{θ} \exp (- x / θ)$ . With a change of variables we can use Gauss Laguerre quadrature

function gauss_laguerre_exponential(N, theta)
    #   E[f(X)] = \int_0^\infty f(x) (1/theta) e^{-x/theta} dx = \int_0^\infty f(theta*y) e^{-y} dy.
    s, w = FastGaussQuadrature.gausslaguerre(N)  # alpha = 0 (default)
    x = theta .* s
    return x, w
end

N = 64
theta = 0.5
x, w = gauss_laguerre_exponential(N, theta)
f(x) = x^2 + 1
@show dot(w, f.(x)), mean(f.(rand(Exponential(theta), 1_000)))

(dot(w, f.(x)), mean(f.(rand(Exponential(theta), 1000)))) = (1.499999999999965, 1.4880275244392855)

(1.499999999999965, 1.4880275244392855)

Similarly, the Gamma distribution with shape parameter $α$ and scale $θ$ has PDF $f (x; α, θ) = \frac{x^{α - 1} e^{- x / θ}}{Γ (α) θ^{α}}$ for $x > 0$ with $Γ (\cdot)$ the Gamma special function.

Using a change of variable and Gauss Laguerre quadrature

function gauss_laguerre_gamma(N, alpha, theta)
  # For Gamma(shape=alpha, scale=theta) with pdf
  #   x^{alpha-1} e^{-x/theta} / (Gamma(alpha) theta^alpha)
  # change variable y = x/theta -> x = theta*y, dx = theta dy
  # E[f(X)] = 1/Gamma(alpha) * ∫_0^∞ f(theta*y) y^{alpha-1} e^{-y} dy
  # FastGaussQuadrature.gausslaguerre(N, a) returns nodes/weights for
  # ∫_0^∞ g(y) y^a e^{-y} dy, so pass a = alpha - 1.

  s, w = FastGaussQuadrature.gausslaguerre(N, alpha - 1)
  x = theta .* s
  w = w ./ SpecialFunctions.gamma(alpha)
  return x, w
end

N = 256
alpha = 7.0
theta = 1.1
x, w = gauss_laguerre_gamma(N, alpha, theta)
f(x) = x^2 + 1
@show dot(w, f.(x)), mean(f.(rand(Gamma(alpha, theta), 100_000)))

(dot(w, f.(x)), mean(f.(rand(Gamma(alpha, theta), 100000)))) = (68.76000000000175, 68.95219068881224)

(68.76000000000175, 68.95219068881224)

7.5. Interpolation #

In economics we often wish to interpolate discrete data (i.e., build continuous functions that join discrete sequences of points).

The package we usually turn to for this purpose is Interpolations.jl.

There are a variety of options, but we will only demonstrate the convenient notations.

7.5.1. Univariate with a Regular Grid#

Let’s start with the univariate case.

We begin by creating some data points, using a sine function

using Interpolations
using Plots

x = -7:7 # x points, coase grid
y = sin.(x) # corresponding y points

xf = -7:0.1:7        # fine grid
plot(xf, sin.(xf), label = "sin function")
scatter!(x, y, label = "sampled data", markersize = 4)

To implement linear and cubic spline interpolation

li = LinearInterpolation(x, y)
li_spline = CubicSplineInterpolation(x, y)

@show li(0.3) # evaluate at a single point

scatter(x, y, label = "sampled data", markersize = 4)
plot!(xf, li.(xf), label = "linear")
plot!(xf, li_spline.(xf), label = "spline")

li(0.3) = 0.25244129544236954

7.5.2. Univariate with Irregular Grid#

In the above, the LinearInterpolation function uses a specialized function for regular grids since x is a Range type.

For an arbitrary, irregular grid

x = log.(range(1, exp(4), length = 10)) .+ 1  # uneven grid
y = log.(x) # corresponding y points

interp = LinearInterpolation(x, y)

xf = log.(range(1,exp(4), length = 100)) .+ 1 # finer grid

plot(xf, interp.(xf), label = "linear")
scatter!(x, y, label = "sampled data", markersize = 4, size = (800, 400))

At this point, Interpolations.jl does not have support for cubic splines with irregular grids, but there are plenty of other packages that do (e.g. Dierckx.jl and GridInterpolations.jl).

7.5.3. Multivariate Interpolation#

Interpolating a regular multivariate function uses the same function

f(x, y) = log(x + y)
xs = 1:0.2:5
ys = 2:0.1:5
A = [f(x, y) for x in xs, y in ys]

# linear interpolation
interp_linear = LinearInterpolation((xs, ys), A)
@show interp_linear(3, 2) # exactly log(3 + 2)
@show interp_linear(3.1, 2.1) # approximately log(3.1 + 2.1)

# cubic spline interpolation
interp_cubic = CubicSplineInterpolation((xs, ys), A)
@show interp_cubic(3, 2) # exactly log(3 + 2)
@show interp_cubic(3.1, 2.1) # approximately log(3.1 + 2.1);

interp_linear(3, 2) = 1.6094379124341003

interp_linear(3.1, 2.1) = 1.6484736801441782
interp_cubic(3, 2) = 1.6094379124341

interp_cubic(3.1, 2.1) =

1.6486586594237707

1.6486586594237707

See Interpolations.jl documentation for more details on options and settings.

7. Quadrature and Interpolation#

7.1. Overview#

7.2. Numerical Integration#