spinterp.RdMonotone interpolation preserves the monotonicity of the data being interpolated, and when the data points are also monotonic, the slopes of the interpolant should also be monotonic.
spinterp(x, y, xp)This implementation follows a cubic version of the method of Delbourgo and Gregory. It yields `shaplier' curves than the Stineman method.
The calculation of the slopes is according to recommended practice:
- monotonic and convex –> harmonic
- monotonic and nonconvex –> geometric
- nonmonotonic and convex –> arithmetic
- nonmonotonic and nonconvex –> circles (Stineman) [not implemented]
The choice of supplementary coefficients r[i] depends on whether
the data are montonic or convex or both:
- monotonic, but not convex
- otherwise
and that can be detected from the data. The choice r[i]=3 for all
i results in the standard cubic Hermitean rational interpolation.
The interpolated values at all the points of xp.
Stan Wagon (2010). Mathematica in Action. Third Edition, Springer-Verlag.
At the moment, the data need to be monotonic and the case of convexity is not considered.
stinepack::stinterp, demography::cm.interp
data1 <- list(x = c(1,2,3,5,6,8,9,11,12,14,15),
y = c(rep(10,6), 10.5,15,50,60,95))
data2 <- list(x = c(0,1,4,6.5,9,10),
y = c(10,4,2,1,3,10))
data3 <- list(x = c(7.99,8.09,8.19,8.7,9.2,10,12,15,20),
y = c(0,0.000027629,0.00437498,0.169183,0.469428,
0.94374,0.998636,0.999919,0.999994))
data4 <- list(x = c(22,22.5,22.6,22.7,22.8,22.9,
23,23.1,23.2,23.3,23.4,23.5,24),
y = c(523,543,550,557,565,575,
590,620,860,915,944,958,986))
data5 <- list(x = c(0,1.1,1.31,2.5,3.9,4.4,5.5,6,8,10.1),
y = c(10.1,8,4.7,4.0,3.48,3.3,5.8,7,7.7,8.6))
data6 <- list(x = c(-0.8, -0.75, -0.3, 0.2, 0.5),
y = c(-0.9, 0.3, 0.4, 0.5, 0.6))
data7 <- list(x = c(-1, -0.96, -0.88, -0.62, 0.13, 1),
y = c(-1, -0.4, 0.3, 0.78, 0.91, 1))
data8 <- list(x = c(-1, -2/3, -1/3, 0.0, 1/3, 2/3, 1),
y = c(-1, -(2/3)^3, -(1/3)^3, -(1/3)^3, (1/3)^3, (1/3)^3, 1))
if (FALSE) { # \dontrun{
opr <- par(mfrow=c(2,2))
# These are well-known test cases:
D <- data1
plot(D, ylim=c(0, 100)); grid()
xp <- seq(1, 15, len=51); yp <- spinterp(D$x, D$y, xp)
lines(spline(D), col="blue")
lines(xp, yp, col="red")
D <- data3
plot(D, ylim=c(0, 1.2)); grid()
xp <- seq(8, 20, len=51); yp <- spinterp(D$x, D$y, xp)
lines(spline(D), col="blue")
lines(xp, yp, col="red")
D <- data4
plot(D); grid()
xp <- seq(22, 24, len=51); yp <- spinterp(D$x, D$y, xp)
lines(spline(D), col="blue")
lines(xp, yp, col="red")
# Fix a horizontal slope at the end points
D <- data8
x <- c(-1.05, D$x, 1.05); y <- c(-1, D$y, 1)
plot(D); grid()
xp <- seq(-1, 1, len=101); yp <- spinterp(x, y, xp)
lines(spline(D, n=101), col="blue")
lines(xp, yp, col="red")
par(opr)} # }