Ajustement de spline cubique multivariée naturelle

17

note: sans réponse correcte après un mois, j'ai republié sur SO

Contexte

J'ai un modèle, , où Y = f ( X )FOui=F(X)

est unematrice n × m d'échantillons à partir de m paramètres et Y est levecteur n × 1 des sorties du modèle.Xn×mmYn×1

est un calcul intensif, donc je voudrais approximer f à l' aide d'une spline cubique multivariée passant par ( X , Y ) points, afin que je puisse évaluer Y à un plus grand nombre de points.ff(X,Y)Y

Question

Existe-t-il une fonction R qui calculera une relation arbitraire entre X et Y?

Plus précisément, je recherche une version multivariée de la splinefunfonction, qui génère une fonction spline pour le cas univarié.

par exemple, voici comment splinefunfonctionne le cas univarié

x <- 1:10
y <- runif(10)
foo <- splinefun(x,y)
foo(1:10) #returns y, as example
all(y == foo(1:10))
## TRUE

Ce que j'ai essayé

J'ai examiné le package mda , et il semble que ce qui suit devrait fonctionner:

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

mais je ne pouvais trouver aucun moyen de mettre en œuvre une spline cubique mars

mise à jour depuis l'offre de la prime, j'ai changé le titre - S'il n'y a pas de fonction R, j'accepterais, par ordre de préférence: une fonction R qui génère une fonction de processus gaussienne, ou une autre fonction d'interpolation multivariée qui passe par les points de conception, de préférence dans R, sinon Matlab.

David LeBauer
la source
essayez la fonction gam (), elle permet n'importe quelle dimension de splines cubiques
user5563

Réponses:

11

Ce document présenté à UseR! 2009 semble résoudre un problème similaire

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

Il suggère le package DiceKriging http://cran.r-project.org/web/packages/DiceKriging/index.html

En particulier, vérifiez les fonctions km et prédisez.

Voici un exemple d'interpolation tridimensionnelle. Il semble simple de généraliser.

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
Tony Ladson
la source
6

Vous avez besoin de plus de données pour un ajustement spline. mgcv est en effet un bon choix. Pour votre demande spécifique, vous devez définir la spline cubique comme fonction de base bs = 'cr' et ne pas la pénaliser avec fx = TRUE. Les deux options sont définies pour un terme fluide défini avec s (). Predict fonctionne comme prévu.

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)
Alex
la source
Merci pour votre aide, mais s'il s'agissait d'une spline cubique, ne devrais-je pas m'attendre predict(foo,x)à revenir y?
David LeBauer
Désolé, je n'ai pas remarqué que vous vouliez une approximation parfaite. Alors, apparemment, mgcv n'est pas d'une grande aide: stop ("La base ne gère que les lissages 1D") (de svn.r-project.org/R-packages/trunk/mgcv/R/smooth.r )
Alex
0

Vous ne donnez aucun détail sur la forme de la fonction F(X); il se peut qu'une fonction constante par morceaux soit une approximation suffisamment bonne, auquel cas vous souhaiterez peut-être ajuster un arbre de régression (avec package rpartpar exemple). Sinon, vous voudrez peut-être regarder le package earth, en plus de ce qui a déjà été suggéré.

F. Tusell
la source
1
la forme de F(X)est au-delà de la portée de ce problème, mais il s'agit d'un modèle de végétation globale dynamique décrit dans l'annexe de Medvigy et al 2009
David LeBauer