评估(即预测)R以外的平滑样条 [英] evaluate (i.e., predict) a smoothing spline outside R
问题描述
我用
library(splines)
Model <- smooth.spline(x, y, df =6)
我想获取拟合的样条,并用外部代码(而不是R中)评估它是否有任意新数据.换句话说,执行predict.smooth.spline
函数的作用.我看了Model
对象:
I would like to take the fitted spline and evaluate it for arbitrary new data in an external code (not in R). In other words, do what the predict.smooth.spline
function does. I had a look at the Model
object:
> 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是有意义的.问题是我不知道如何使用Model$fit$knot
和Model$fit$coef
在R之外进行预测.我试图看一下predict.smooth.spline
,但令人惊讶的是,我得到了
I think the Model$fit$knot
and Model$fit$coef
vectors contain the full description of the fit. Note that the knots are 20, while x
and y
have 14 elements each: I always thought a smoothing spline would have as many knots as fitting points. However, since the first three and the last three knots are identical, 20-6 = 14 which makes sense. The problem is that I don't know how to use Model$fit$knot
and Model$fit$coef
to make predictions outside R. I tried to have a look at predict.smooth.spline
, but surprisingly that's what I get
> library(splines)
> predict.smooth.spline
Error: object 'predict.smooth.spline' not found
由于显然有些用户误解了这个问题,我知道如何在R中使用predict
来获取我的平滑样条曲线的新值.问题是我想在外部代码中做出这些预测.因此,我想看一下函数predict.smooth.spline
的代码,以便尝试在R之外重现该算法.通常,在R中,您只需输入函数名即可读取函数代码(不带参数,不带参数).括号).但是,当我尝试使用predict.smooth.spline
执行此操作时,出现了以上错误.
since apparently some users misunderstood the question, I know how to use predict
in R, to get new values of my smoothing spline. The problem is that I want to make those predictions in an external code. Thus I wanted to have a look at the code for the function predict.smooth.spline
, so that I could try to reproduce the algorithm outside R. Usually in R you can read the code of a function just by entering its name (without arguments and without parentheses) at the R prompt. But when I try to do that with predict.smooth.spline
, I get the above error.
感谢@ r2evans的大力帮助,我找到了smooth.spline
的predict
方法的源代码.我(认为我)了解其中的大部分内容:
thanks to the great help from @r2evans, I found the source for the predict
method of smooth.spline
. I (think I) understand most of it:
> 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)
}
但是,我有两个困难:
the
.Fortran()
function calls a Fortran subroutinebvalus
whose source is quite simple. However,bvalus
in turn callsbvalue
which is much more complicated, and callsinterv
whose source I cannot find. Bad news:bvalue
is way too complicated for me to understand (I'm definitely not a Fortran expert). Good news: the external code which must reproducepredict.smooth.spline.fit
is also a Fortran code. If worse comes to worst, I could just ask my coworker to include the source frombvalus
andbvalue
in his code. However, even in this admittedly not so nice scenario, I would still miss the source code forinterv
(I hope it doesn't call something else!!!).
我不明白这里正在做什么(请注意,我只对deriv == 0
案例感兴趣):
I don't understand what it's being done here (note I'm only interested in the deriv == 0
case):
k
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)
}
某种递归代码?这里有帮助吗?
Some sort of recursive code? Any help here?
推荐答案
smooth.spline
不在splines
软件包中,而是在stats
中.此外,它不会导出,因此您必须使用三冒号方法才能看到它:stats:::predict.smooth.spline
.然后将您指向predict.smooth.spline.fit
,可以通过类似的方式找到它. (由于它可选地使用.Fortran()
,因此您可能必须推断正在发生的事情……除非您深入研究
smooth.spline
is not in the splines
package, it's in stats
. Additionally, it's not exported, so you have to use the triple-colon method to see it: stats:::predict.smooth.spline
. It then points you to predict.smooth.spline.fit
, which can be found in a similar manner. (Since it optionally uses .Fortran()
, you may have to infer what is going on ... unless you dive into the source.)
这篇关于评估(即预测)R以外的平滑样条的文章就介绍到这了,希望我们推荐的答案对大家有所帮助,也希望大家多多支持IT屋!