拟合多元、自然三次样条

机器算法验证 r 多元分析 样条 插值 高斯过程
2022-01-25 11:53:54

注意:一个月后没有正确答案,我已重新发布到SO

背景

我有一个模型,f, 在哪里Y=f(X)

X是一个n×m样本矩阵m参数和Y是个n×1模型输出向量。

f是计算密集型的,所以我想近似f使用多元三次样条通过(X,Y)点,以便我可以评估Y在更多的点上。

问题

是否有一个 R 函数可以计算 X 和 Y 之间的任意关系?

具体来说,我正在寻找该splinefun函数的多变量版本,它为单变量情况生成样条函数。

例如,这就是splinefun单变量案例的工作方式

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

我试过的

我已经查看了mda包,似乎以下内容应该有效:

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

但我找不到任何方法来实现三次样条mars

自从提供赏金以来更新,我更改了标题 - 如果没有 R 函数,我会按优先顺序接受:输出高斯过程函数的 R 函数,或通过设计点的另一个多元插值函数,最好在 R 中,否则在 Matlab 中。

3个回答

这篇论文发表在 UserR!2009 似乎解决了类似的问题

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

它建议 DiceKriging 包http://cran.r-project.org/web/packages/DiceKriging/index.html

特别是检查函数 km 和 predict。

这是一个三维插值的例子。它看起来很容易概括。

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

您需要更多数据来进行样条拟合。mgcv 确实是一个不错的选择。对于您的特定要求,您需要将三次样条设置为基函数 bs='cr' 并且不要用 fx=TRUE 对其进行惩罚。这两个选项都设置为使用 s() 设置的平滑项。预测按预期工作。

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)

您没有提供有关功能形式的详细信息f(X); 可能分段常数函数是一个足够好的近似值,在这种情况下,您可能想要拟合回归树(rpart例如使用包)。earth否则,除了已经建议的内容之外,您可能还想查看 package 。