# Fitting multivariate, natural cubic spline

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

Background

I have a model, $f$, where $Y=f(\textbf{X})$

$\textbf{X}$ is an $n \times m$ matrix of samples from $m$ parameters and $Y$ is the $n \times 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.

Question

Is 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 tried

I 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

update 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.

-

You need more data for a spline fit. mgcv indeed is a good choice. For your specific request you need to set the cubic spline as the basis function bs='cr' and also not have it penalized with fx=TRUE. Both options are set for a smooth term that is set with s(). Predict works as expected.

library(mgcv)
x <- data.frame(a = 1:100, b = 1:100/2, c = 1:100*2)
y <- runif(100)
foo <- gam(y~a+b+s(c,bs="cr",fx=TRUE),data=x)
plot(foo)
predict(foo,x)

-
 Thank you for your help, but if this were a cubic spline, shouldn't I expect predict(foo,x) to return y? – David Jul 28 '11 at 21:57 Sorry, didn't notice that you want a perfect approximation. Then apparently mgcv is not of much help: stop("Basis only handles 1D smooths") (from svn.r-project.org/R-packages/trunk/mgcv/R/smooth.r) – Alex Jul 29 '11 at 0:38

try gam() function, it allows any dimension of cubic splines

-
 Maybe you could elaborate a bit? – mbq♦ Jul 27 '11 at 8:26 @mbq I don't remember all details. There are two packages - gam and mgcv, both having the gam() function. I worked with mgcv package and know that it allows different splines - thin-plate, cubic and more. Here is its help: link, see pp.28-30 and 46 – user5563 Jul 27 '11 at 10:37 @user5563 the gam::gam function does not support cubic spline, and I am not able to get mgcv::gam to do so. Could you please provide an example? – David Aug 4 '11 at 6:01 I've not had the gam() from the gam package work for me, while the one from the mgcv package works great. It might only be me, but I'd recommend mgcv myself. – Wayne Aug 22 '11 at 17:28

You give no details as to the form of the function $f(X)$; it might be that a piecewise constant function is a sufficiently good approximation, in which case you might want to fit a regression tree (with package rpart for instance). Otherwise, you might want to look at package earth, in addition to what has been suggested already.

-
 the form of $f(X)$ is beyond the scope of this problem, but it is a dynamic global vegetation model described in the appendix of Medvigy et al 2009 – David Aug 23 '11 at 4:18