0

我将平滑样条曲线拟合到 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$knotModel$fit$coef向量包含拟合的完整描述。请注意,节点有 20 个,而x每个节点y有 14 个元素:我一直认为平滑样条曲线的节点数与拟合点一样多。但是,由于前三个和后三个结是相同的,所以 20-6 = 14 这是有道理的。问题是我不知道如何在 R 之外使用Model$fit$knotModel$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 的大力帮助,我找到predictsmooth.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)
}

但是,我有两个困难:

  1. .Fortran()函数调用一个源代码非常简单bvalus的Fortran 子例程。但是,依次调用更复杂的调用,以及我找不到源的调用。坏消息:对我来说太复杂了(我绝对不是 Fortran 专家)。好消息:必须重现的外部代码也是 Fortran 代码。如果情况变得更糟,我可以要求我的同事在他的代码中包含源代码。然而,即使在这个公认的不太好的场景中,我仍然会错过源代码(我希望它不会调用别的东西!!!)。bvalusbvalueintervbvaluepredict.smooth.spline.fitbvalusbvalueinterv

  2. 我不明白这里做了什么(注意我只对这个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)
        }

某种递归代码?这里有什么帮助吗?

4

2 回答 2

2

smooth.spline不在splines包装内,在stats. 另外,它没有导出,所以你必须使用三冒号方法才能看到它:stats:::predict.smooth.spline. 然后它会将您指向predict.smooth.spline.fit,可以以类似的方式找到它。(由于它可以选择使用.Fortran(),您可能必须推断发生了什么......除非您深入了解源代码。)

于 2016-02-27T02:55:59.033 回答
1

将平滑样条导出为分段多项式是在 R 之外重建样条的一种方法。我的包SplinesUtilshttps ://github.com/ZheyuanLi/SplinesUtils可以做到这一点。你可以通过

devtools::install_github("ZheyuanLi/SplinesUtils")

这里要使用的函数是SmoothSplinesAsPiecePoly. 我只是在其文档中复制此功能的示例。

library(SplinesUtils)

## a toy dataset
set.seed(0)
x <- 1:100 + runif(100, -0.1, 0.1)
y <- poly(x, 9) %*% rnorm(9)
y <- y + rnorm(length(y), 0, 0.2 * sd(y))

## fit a smoothing spline
sm <- smooth.spline(x, y)

## coerce "smooth.spline" object to "PiecePoly" object
oo <- SmoothSplineAsPiecePoly(sm)

## print the "PiecePoly"
oo
#61 piecewise polynomials of degree 3 are constructed!
#Use 'summary' to export all of them.
#The first 6 are printed below.
#-0.626 - 0.17 * (x - 1.08) - 0 * (x - 1.08) ^ 2 - 0.0094 * (x - 1.08) ^ 3
#-0.768 - 0.148 * (x - 1.95) - 0.0246 * (x - 1.95) ^ 2 - 0.00569 * (x - 1.95) ^ 3
#-0.919 + 0.0259 * (x - 4.01) + 0.0598 * (x - 4.01) ^ 2 + 0.0086 * (x - 4.01) ^ 3
#-0.834 + 0.124 * (x - 5.08) + 0.0323 * (x - 5.08) ^ 2 + 0.00466 * (x - 5.08) ^ 3
#-0.494 + 0.197 * (x - 7.08) + 0.00433 * (x - 7.08) ^ 2 + 0.0027 * (x - 7.08) ^ 3
#-0.113 + 0.183 * (x - 9.03) + 0.0115 * (x - 9.03) ^ 2 + 0.00377 * (x - 9.03) ^ 3

样条的结是

oo$knots
# [1]   1.079339   1.953102   4.014571   5.081642   7.079678   9.032160
# [7]  10.025823  11.941195  12.935311  14.976821  16.999540  18.043524
#[13]  19.976007  22.086941  22.942429  24.925111  25.953444  27.902678
#[19]  30.073938  30.968070  33.019913  34.937244  36.065475  38.058848
#[25]  38.921589  40.982255  43.029412  44.056587  46.005944  47.904666
#[31]  48.995446  51.038546  51.995524  53.987619  55.914136  56.919893
#[37]  59.003727  60.981366  62.082575  63.991813  64.966479  66.951603
#[43]  69.053262  69.916849  71.967815  73.969337  74.966755  77.078440
#[49]  78.072868  80.055464  81.986932  83.042503  84.965070  86.940538
#[55]  88.042224  89.949098  90.928661  92.911787  95.075254  96.055783
#[61]  97.991055 100.062174

分段多项式系数是

CoefMat <- oo$PiecePoly$coef
## for the first 5 pieces
CoefMat[, 1:5]
#             [,1]         [,2]         [,3]        [,4]         [,5]
#[1,] -0.626225456 -0.768321912 -0.919380838 -0.83408278 -0.494257767
#[2,] -0.169805245 -0.148267472  0.025888868  0.12418698  0.197353171
#[3,]  0.000000000  0.024649465  0.059832206  0.03228737  0.004331680
#[4,]  0.009403574  0.005688943 -0.008604501 -0.00466386 -0.002702257

该软件包还具有其他功能。有关示例,请参阅通过 R 函数“smooth.spline”识别拟合平滑样条的所有局部极值。

于 2018-10-08T00:23:58.277 回答