Interpolation
Interpolation is the process of constructing a function that takes on specified values at specified points. In engineering and science, one often has a number of data points, obtained by sampling or experimentation, which represent the values of a function for a limited number of values of the independent variable. It is often required to interpolate (i.e. estimate) the value of that function for an intermediate value of the independent variable.
Smile package smile.interpolation provides a variety of
algorithms on 1d and 2d data. These algorithms implements the interface
Interpolation (or Interpolation2D for 2d data),
of which the method interpolate takes a value and return
an interploated value.
Piecewise Linear Interpolation
LinearInterpolation is quick and easy, but it is not very precise. Another disadvantage
is that the interpolant is not differentiable at the control points.
val x = Array(0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0)
val y = Array(0.0, 0.8415, 0.9093, 0.1411, -0.7568, -0.9589, -0.2794)
val linear = new LinearInterpolation(x, y)
val data = (0 to 60).map { i =>
val x = i * 0.1
val y = linear.interpolate(x)
Array(x, y)
}.toArray
val canvas = line(data)
show(canvas)
double[] x = {0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0}
double[] y = {0.0, 0.8415, 0.9093, 0.1411, -0.7568, -0.9589, -0.2794}
var linear = new LinearInterpolation(x, y)
double[][] data = new double[61][2]
for (int i = 0; i < data.length; i++) {
data[i][0] = i * 0.1;
data[i][1] = linear.interpolate(data[i][0]);
}
var canvas = LinePlot.plot(data, Line.Style.SOLID)
canvas.window()
BilinearInterpolation is an extension of linear interpolation for
interpolating functions of two variables on a
regular grid. The key idea is to perform linear interpolation first in one
direction, and then again in the other direction.
Cubic Spline Interpolation
Spline interpolation uses low-degree polynomials in each of the intervals, and chooses the polynomial pieces such that they fit smoothly together. The resulting function is called a spline. The natural cubic spline is piecewise cubic and twice continuously differentiable. Furthermore, its second derivative is zero at the end points.
Like polynomial interpolation, spline interpolation incurs a smaller error than linear interpolation and the interpolant is smoother. However, the interpolant is easier to evaluate than the high-degree polynomials used in polynomial interpolation. It also does not suffer from Runge's phenomenon.
val cubic = new CubicSplineInterpolation1D(x, y)
val data = (0 to 60).map { i =>
val x = i * 0.1
val y = cubic.interpolate(x)
Array(x, y)
}.toArray
canvas.line("Cubic", data, RED)
show(canvas)
var cubic = new CubicSplineInterpolation1D(x, y)
double[][] data = new double[61][2]
for (int i = 0; i < data.length; i++) {
data[i][0] = i * 0.1;
data[i][1] = cubic.interpolate(data[i][0]);
}
canvas.line("Cubic", data, Color.RED)
canvas.window()
For 2d grid, we provide two implementations: CubicSplineInterpolation2D and
BicubicInterpolation. CubicSplineInterpolation2D is
similar to one-dimensional splines as it guarantees the
continuity of the first and second function derivatives. In contrast,
BicubicInterpolation guarantees continuity of only
gradient and cross-derivative. Second derivatives
could be discontinuous.
In image processing, bicubic interpolation is often chosen over bilinear interpolation or nearest neighbor in image resampling, when speed is not an issue. Images resampled with bicubic interpolation are smoother and have fewer interpolation artifacts.
val x1 = Array(1950.0, 1960, 1970, 1980, 1990)
val x2 = Array(10.0, 20, 30)
val y = Array(
Array(150.697, 199.592, 187.625),
Array(179.323, 195.072, 250.287),
Array(203.212, 179.092, 322.767),
Array(226.505, 153.706, 426.730),
Array(249.633, 120.281, 598.243)
)
val canvas = heatmap(y, Palette.jet(256))
canvas.setTitle("Original")
show(canvas)
val cubic = new CubicSplineInterpolation2D(x1, x2, y)
val data = Array.ofDim[Double](101, 101)
for (i <- 0 to 100; j <- 0 to 100)
data(i)(j) = cubic.interpolate(1950 + i*0.4, 10 + j*0.2)
val cubicPlot = heatmap(data, Palette.jet(256))
cubicPlot.setTitle("Cubic")
show(cubicPlot)
val bicubic = new BicubicInterpolation(x1, x2, y)
for (i <- 0 to 100; j <- 0 to 100)
data(i)(j) = bicubic.interpolate(1950 + i*0.4, 10 + j*0.2)
val bicubicPlot = heatmap(data, Palette.jet(256))
bicubicPlot.setTitle("Bicubic")
show(bicubicPlot)
double[] x1 = {1950.0, 1960, 1970, 1980, 1990}
double[] x2 = {10.0, 20, 30}
double[][] y = {
{150.697, 199.592, 187.625},
{179.323, 195.072, 250.287},
{203.212, 179.092, 322.767},
{226.505, 153.706, 426.730},
{249.633, 120.281, 598.243}
}
var canvas = Heatmap.plot(y, Palette.jet(256))
canvas.setTitle("Original")
canvas.window()
var cubic = new CubicSplineInterpolation2D(x1, x2, y)
double[][] data = new double[101][101]
for (int i = 0; i <= 100; i++) {
for (int j = 0; j <= 100; j++) {
data[i][j] = cubic.interpolate(1950 + i*0.4, 10 + j*0.2);
}
}
var cubicPlot = Heatmap.plot(data, Palette.jet(256))
cubicPlot.setTitle("Cubic")
cubicPlot.window()
val bicubic = new BicubicInterpolation(x1, x2, y)
for (int i = 0; i <= 100; i++) {
for (int j = 0; j <= 100; j++) {
data[i][j] = cubic.interpolate(1950 + i*0.4, 10 + j*0.2);
}
}
var bicubicPlot = Heatmap.plot(data, Palette.jet(256))
bicubicPlot.setTitle("Bicubic")
bicubicPlot.window()
Kriging Interpolation
Kriging interpolation can be used for the data points irregularly distributed in space. Kriging belongs to the family of linear least squares estimation algorithms, also known as Gauss-Markov estimation or Gaussian process regression. We implement ordinary kriging for interpolation with power variogram.
val kriging = new KrigingInterpolation1D(x, y)
val data = Array.ofDim[Double](61, 2)
for (i <- 0 to 60) {
data(i)(0) = i * 0.1
data(i)(1) = kriging.interpolate(data(i)(0))
}
canvas.line("Kriging", data, CYAN)
show(canvas)
var kriging = new KrigingInterpolation1D(x, y)
double[][] data = new double[61][2]
for (int i = 0; i < data.length; i++) {
data[i][0] = i * 0.1;
data[i][1] = kriging.interpolate(data[i][0]);
}
canvas.line("Kriging", data, Color.CYAN)
canvas.window()
RBF Interpolation
Radial basis function interpolation is a popular method for the data points are
irregularly distributed in space. In its basic form, radial basis function
interpolation is in the form
y(x) = Σ wi φ(||x-ci||),
where the approximating function y(x) is represented as a sum of N radial
basis functions φ, each associated with a different center ci, and weighted
by an appropriate coefficient wi. For distance, one usually chooses
euclidean distance. The weights wi can
be estimated using the matrix methods of linear least squares, because
the approximating function is linear in the weights.
The points ci often called the centers or collocation points of the RBF interpolant. Note also that the centers ci can be located at arbitrary points in the domain, and do not require a grid. For certain RBF exponential convergence has been shown. Radial basis functions were successfully applied to problems as diverse as computer graphics, neural networks, for the solution of differential equations via collocation methods and many other problems.
Other popular choices for φ comprise the Gaussian function and the so called thin plate splines. Thin plate splines result from the solution of a variational problem. The advantage of the thin plate splines is that their conditioning is invariant under scalings. Gaussians, multi-quadrics and inverse multi-quadrics are infinitely smooth and and involve a scale or shape parameter, r0 > 0. Decreasing r0 tends to flatten the basis function. For a given function, the quality of approximation may strongly depend on this parameter. In particular, increasing r0 has the effect of better conditioning (the separation distance of the scaled points increases).
A variant on RBF interpolation is normalized radial basis function (NRBF) interpolation, in which we require the sum of the basis functions to be unity. NRBF arises more naturally from a Bayesian statistical perspective. However, there is no evidence that either the NRBF method is consistently superior to the RBF method, or vice versa.
val rbf = new RBFInterpolation1D(x, y, new smile.math.rbf.GaussianRadialBasis())
val data = Array.ofDim[Double](61, 2)
for (i <- 0 to 60) {
data(i)(0) = i * 0.1
data(i)(1) = rbf.interpolate(data(i)(0))
}
canvas.line("RBF", data, GREEN)
show(canvas)
var rbf = new RBFInterpolation1D(x, y, new smile.math.rbf.GaussianRadialBasis())
double[][] data = new double[61][2]
for (int i = 0; i < data.length; i++) {
data[i][0] = i * 0.1;
data[i][1] = rbf.interpolate(data[i][0]);
}
canvas.line("RBF", data, Color.GREEN)
canvas.window()
Shepard Interpolation
Shepard interpolation is a special case of normalized radial basis function interpolation if the function φ(r) goes to infinity as r → 0, and is finite for r > 0. In this case, the weights wi are just equal to the respective function values yi. So we need not solve linear equations and thus it works for very large N. An example of such φ is φ(r) = r-p with (typically) 1 < p ≤ 3.
val shepard = new ShepardInterpolation1D(x, y, 3)
val data = Array.ofDim[Double](61, 2)
for (i <- 0 to 60) {
data(i)(0) = i * 0.1
data(i)(1) = shepard.interpolate(data(i)(0))
}
canvas.line("Shepard", data, PINK)
show(canvas)
var shepard = new ShepardInterpolation1D(x, y, 3)
double[][] data = new double[61][2]
for (int i = 0; i < data.length; i++) {
data[i][0] = i * 0.1;
data[i][1] = shepard.interpolate(data[i][0]);
}
canvas.line("Shepard", data, Color.PINK)
canvas.window()
Shepard interpolation is rarely as accurate as the well-tuned application of other radial basis functions. However, it is simple, fast, and often sufficient for quick and dirty applications.
Laplace Interpolation
Laplace interpolation can restore missing or unmeasured values on a 2-dimensional evenly spaced regular grid. In some sense, Laplace interpolation produces the smoothest possible interpolant, which are obtained by solving a very sparse linear equations with biconjugate gradient method.
val zz = Array.ofDim[Double](101, 101)
val ww = Array.ofDim[Double](101, 101)
for (i <- 0 to 100; j <- 0 to 100) {
zz(i)(j) = if (java.lang.Math.random() < 0.2) Double.NaN else data(i)(j)
ww(i)(j) = zz(i)(j)
}
val canvas = heatmap(ww, Palette.jet(256))
canvas.setTitle("Missing Values")
show(canvas)
LaplaceInterpolation.interpolate(zz)
val canvas = heatmap(zz, Palette.jet(256))
canvas.setTitle("Laplace")
show(canvas)
double[][] zz = new double[101][101]
double[][] ww = new double[101][101]
for (int i = 0; i <= 100; i++) {
for (int j = 0; j <= 100; j++) {
zz[i][j] = Math.random() < 0.2 ? Double.NaN : data[i][j];
ww[i][j] = zz[i][j];
}
}
var canvas = Heatmap.plot(ww, Palette.jet(256))
canvas.setTitle("Missing Values")
canvas.window()
LaplaceInterpolation.interpolate(zz)
var canvas = Heatmap.plot(zz, Palette.jet(256))
canvas.setTitle("Laplace")
canvas.window()
