note:with no correct answers after a month, I have reposted to SO

BackgroundI have a model, f, where Y=f(X)

X is an n×m matrix of samples from m parameters and Y is the n×1 vector of model outputs.

f is computationally intensive, so I would like to approximate f using a multivariate cubic spline through (X,Y) points, so that I can evaluate Y at a larger number of points.

QuestionIs there an R function that will calculate an arbitrary relationship between X and Y?

Specifically, I am looking for a multivariate version of the

`splinefun`

function, which generates a spline function for the univariate case.e.g. this is how

`splinefun`

works for the univariate case`x <- 1:10 y <- runif(10) foo <- splinefun(x,y) foo(1:10) #returns y, as example all(y == foo(1:10)) ## TRUE`

What I have triedI have reviewed the mda package, and it seems that the following should work:

`library(mda) x <- data.frame(a = 1:10, b = 1:10/2, c = 1:10*2) y <- runif(10) foo <- mars(x,y) predict(foo, x) #all the same value all(y == predict(foo,x)) ## FALSE`

but I could not find any way to implement a cubic-spline in

`mars`

since offering the bounty, I changed the title – If there is no R function, I would accept, in order of preference: an R function that outputs a gaussian process function, or another multivariate interpolating function that passes through the design points, preferably in R, else Matlab.update

**Answer**

This paper presented at UseR! 2009 seems to address a similar problem

http://www.r-project.org/conferences/useR-2009/slides/Roustant+Ginsbourger+Deville.pdf

It suggests the DiceKriging package http://cran.r-project.org/web/packages/DiceKriging/index.html

In particular, check the functions km and predict.

Here is a an example of three dimensional interpolation. It looks to be straightforward to generalise.

```
x <- c(0, 0.4, 0.6, 0.8, 1)
y <- c(0, 0.2, 0.3, 0.4, 0.5)
z <- c(0, 0.3, 0.4, 0.6, 0.8)
model <- function(param){
2*param[1] + 3*param[2] +4*param[3]
}
model.in <- expand.grid(x,y,z)
names(model.in) <- c('x','y','z')
model.out <- apply(model.in, 1, model)
# fit a kriging model
m.1 <- km(design=model.in, response=model.out, covtype="matern5_2")
# estimate a response
interp <- predict(m.1, newdata=data.frame(x=0.5, y=0.5, z=0.5), type="UK", se.compute=FALSE)
# check against model output
interp$mean
# [1] 4.498902
model(c(0.5,0.5,0.5))
# [1] 4.5
# check we get back what we put in
interp <- predict(m.1, newdata=model.in, type="UK", se.compute=FALSE)
all.equal(model.out, interp$mean)
# TRUE
```

**Attribution***Source : Link , Question Author : David LeBauer , Answer Author : Tony Ladson*