我将平滑样条曲线拟合到 R 中的数据
library(splines)
Model <- smooth.spline(x, y, df =6)
我想采用拟合的样条曲线并评估它以获取外部代码(而不是 R)中的任意新数据。换句话说,做predict.smooth.spline
函数所做的事情。我看了看Model
对象:
> str(Total_work_model)
List of 15
$ x : num [1:14] 0.0127 0.0186 0.0275 0.0343 0.0455 ...
$ y : num [1:14] 3174 3049 2887 2862 2975 ...
$ w : num [1:14] 1 1 1 1 1 1 1 1 1 1 ...
$ yin : num [1:14] 3173 3075 2857 2844 2984 ...
$ data :List of 3
..$ x: num [1:14] 0.0343 0.0455 0.0576 0.0697 0.0798 ...
..$ y: num [1:14] 2844 2984 3048 2805 2490 ...
..$ w: num [1:14] 1 1 1 1 1 1 1 1 1 1 ...
$ lev : num [1:14] 0.819 0.515 0.542 0.568 0.683 ...
$ cv.crit : num 6494075
$ pen.crit: num 3260
$ crit : num 3
$ df : num 8
$ spar : num 0.353
$ lambda : num 8.26e-05
$ iparms : Named int [1:3] 3 0 10
..- attr(*, "names")= chr [1:3] "icrit" "ispar" "iter"
$ fit :List of 5
..$ knot : num [1:20] 0 0 0 0 0.056 ...
..$ nk : int 16
..$ min : num 0.0127
..$ range: num 0.104
..$ coef : num [1:16] 3174 3132 3027 2871 2842 ...
..- attr(*, "class")= chr "smooth.spline.fit"
$ call : language smooth.spline(x = Phi00, y = Total, df = 8)
- attr(*, "class")= chr "smooth.spline"
我认为Model$fit$knot
和Model$fit$coef
向量包含拟合的完整描述。请注意,节点有 20 个,而x
每个节点y
有 14 个元素:我一直认为平滑样条曲线的节点数与拟合点一样多。但是,由于前三个和后三个结是相同的,所以 20-6 = 14 这是有道理的。问题是我不知道如何在 R 之外使用Model$fit$knot
和Model$fit$coef
做出预测。我试图看看predict.smooth.spline
,但令人惊讶的是我得到了
> library(splines)
> predict.smooth.spline
Error: object 'predict.smooth.spline' not found
编辑:由于显然有些用户误解了这个问题,我知道如何predict
在 R 中使用,以获得我的平滑样条曲线的新值。问题是我想在外部代码中做出这些预测。因此我想看看函数的代码predict.smooth.spline
,这样我就可以尝试在 R 之外重现算法。通常在 R 中,您只需在R 提示符。但是当我尝试这样做时predict.smooth.spline
,我得到了上述错误。
EDIT2:感谢@r2evans 的大力帮助,我找到predict
了smooth.spline
. 我(想我)了解其中的大部分内容:
> stats:::predict.smooth.spline.fit
function (object, x, deriv = 0, ...)
{
if (missing(x))
x <- seq.int(from = object$min, to = object$min + object$range,
length.out = length(object$coef) - 4L)
xs <- (x - object$min)/object$range
extrap.left <- xs < 0
extrap.right <- xs > 1
interp <- !(extrap <- extrap.left | extrap.right)
n <- sum(interp)
y <- xs
if (any(interp))
y[interp] <- .Fortran(C_bvalus, n = as.integer(n), knot = as.double(object$knot),
coef = as.double(object$coef), nk = as.integer(object$nk),
x = as.double(xs[interp]), s = double(n), order = as.integer(deriv))$s
if (any(extrap)) {
xrange <- c(object$min, object$min + object$range)
if (deriv == 0) {
end.object <- Recall(object, xrange)$y
end.slopes <- Recall(object, xrange, 1)$y * object$range
if (any(extrap.left))
y[extrap.left] <- end.object[1L] + end.slopes[1L] *
(xs[extrap.left] - 0)
if (any(extrap.right))
y[extrap.right] <- end.object[2L] + end.slopes[2L] *
(xs[extrap.right] - 1)
}
else if (deriv == 1) {
end.slopes <- Recall(object, xrange, 1)$y * object$range
y[extrap.left] <- end.slopes[1L]
y[extrap.right] <- end.slopes[2L]
}
else y[extrap] <- 0
}
if (deriv > 0)
y <- y/(object$range^deriv)
list(x = x, y = y)
}
但是,我有两个困难:
该
.Fortran()
函数调用一个源代码非常简单bvalus
的Fortran 子例程。但是,依次调用更复杂的调用,以及我找不到源的调用。坏消息:对我来说太复杂了(我绝对不是 Fortran 专家)。好消息:必须重现的外部代码也是 Fortran 代码。如果情况变得更糟,我可以要求我的同事在他的代码中包含源代码。然而,即使在这个公认的不太好的场景中,我仍然会错过源代码(我希望它不会调用别的东西!!!)。bvalus
bvalue
interv
bvalue
predict.smooth.spline.fit
bvalus
bvalue
interv
我不明白这里做了什么(注意我只对这个
deriv == 0
案例感兴趣):
ķ
if (any(extrap)) {
xrange <- c(object$min, object$min + object$range)
if (deriv == 0) {
end.object <- Recall(object, xrange)$y
end.slopes <- Recall(object, xrange, 1)$y * object$range
if (any(extrap.left))
y[extrap.left] <- end.object[1L] + end.slopes[1L] *
(xs[extrap.left] - 0)
if (any(extrap.right))
y[extrap.right] <- end.object[2L] + end.slopes[2L] *
(xs[extrap.right] - 1)
}
某种递归代码?这里有什么帮助吗?