smint (version 0.4.0)

qsheppInt: Shepard's (modified) quadratic interpolation method

Description

Shepard's (modified) quadratic interpolation method.

Usage

qsheppInt(X, y, XNew = NULL, nQ = 15L, nW = 33L, nR,
          checkX = TRUE)
qsheppInt2d(X, y, XNew = NULL, nQ = 13L, nW = 19L, nR,
            deriv = 0L, checkX = TRUE)
qsheppInt3d(X, y, XNew = NULL, nQ = 17L, nW = 32L, nR,
            deriv = 0L, checkX = TRUE)

Arguments

X

Design matrix. A matrix with $n$ rows and $d$ columns where $n$ is the number of nodes and $d$ is the dimension.

y

Numeric vector of length $n$ containing the function values at design points.

XNew

Numeric matrix of new design points where interpolation is computed. It must have $d$ columns.

nQ

Number of points used in the local polynomial fit. This is the parameter NQ of the Fortran routine qshepmd and the parameter $N_p$ of the referenced article.

nW

Number of nodes within (and defining) the radii of influence, which enter into the weights. This is parameter NW of the Fortran routine qshepmd and $N_w$ of the reference article.

nR

Number of divisions in each dimension for the cell grid defined in the subroutine storem. A hyperbox containing the nodes is partitioned into cells in order to increase search efficiency. The recommended value $(n/3)^(1/d)$ is used as defau

checkX

If TRUE, the dimensions and colnames of X and XNew will be checked.

deriv

Logical or integer value in c(0, 1). When the (coerced) integer value is 1, the derivatives are computed. This is is only possible for the dimensions 2 and 3, and only when the specific interpolation fun

Value

• A list with several objects related to the computation method. The vector of interpolated value is in the list element named yNew.

Details

Shepard's modified interpolation method works for scattered data in arbitrary dimension. It relies on a local polynomial fit and the quadratic version uses a polynomial of degree 2 (a quadratic form) as local approximation of the function.

References

W.I. Thacker, J. Zhang, L.T. Watson, J.B. Birch, M.A. Iyer and M.W. Berry (2010). Algorithm 905: SHEPPACK: Modified Shepard Algorithm for Interpolation of Scattered Multivariate Data ACM Trans. on Math. Software (TOMS) Vol. 37, n. 3. http://dl.acm.org/citation.cfm?id=1824812{link}

M.W. Berry and K.S. Minser (1999). Algorithm 798: High-dimensional interpolation using the modified Shepard method. ACM Trans. Math. Software (TOMS) Vol. 25, n. 3, pp. 353-366. http://dl.acm.org/citation.cfm?id=326147.326154{link}

Examples

n <- 1500; nNew <- 100; d <- 4
fTest <- function(x)((x[1] + 2 * x[2]   + 3 * x[3] + 4 * x[4]) / 12)^2
set.seed(12345)
X <- matrix(runif(n*d), nrow = n, ncol = d)
y <- apply(X, 1, FUN = fTest)
XNew <- matrix(runif(nNew * d), nrow = nNew, ncol = d)
yNew <- apply(XNew, 1, FUN = fTest)
system.time(res <- qsheppInt(X = X, XNew = XNew, y = y, nQ = 40,
                             checkX = FALSE))
## check errors
max(abs(res$yNew - yNew))

##=========================================================================
## Use SHEPPACK test functions see Thacker et al. section 7 'PERFORMANCE'
##=========================================================================
set.seed(1234)
   d <- 3
   k <- 0:4; n0 <- 100 * 2^k; n1 <- 4
   GD <- Grid(nlevels = rep(n1, d))
   XNew <- as.matrix(GD)
   RMSE <- array(NA, dim = c(5, length(k)),
                 dimnames = list(fun = 1:5, k = k))
   for (iFun in 1:5) {
      yNew <- apply(XNew, 1, ShepFuns[[iFun]])
      for (iN in 1:length(n0)) {
         X <- matrix(runif(n0[iN] * d), ncol = d)
         y <- apply(X, 1, ShepFuns[[iFun]])
         res <- qsheppInt(X = X, XNew = XNew, y = y, nQ = 40, checkX = FALSE)
         RMSE[iFun, iN] <- mean((res$yNew - yNew)^2)
      }
   }
   cols <- c("black", "SteelBlue2", "orangered", "SpringGreen3", "purple")
   pchs <-  c(16, 21, 22, 23, 24)
   matplot(k, t(RMSE), type = "o", lwd = 2, lty = 1,
           col = cols, xaxt = "n", pch = pchs, cex = 1.4,
           bg = "white",
           main = sprintf("dim = %d SHEPPACK test functions", d),
           xlab = "number of nodes", ylab = "RMSE")
   axis(side = 1, at = k, labels = n0)
   legend("topright", legend = paste("shepFun", 1:5),
           col = cols, pch = pchs, lwd = 2, pt.lwd = 2, pt.bg = "white")