Ajusté una spline de suavizado a los datos en R con

library(splines)
Model <- smooth.spline(x, y, df =6) 

Me gustaría tomar la spline ajustada y evaluarla en busca de nuevos datos arbitrarios en un código externo (no en R). En otras palabras, haz lo quepredict.smooth.spline función hace. Eché un vistazo a laModel objeto:

> 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"

Pienso que elModel$fit$knot yModel$fit$coef Los vectores contienen la descripción completa del ajuste. Tenga en cuenta que los nudos son 20, mientras quex yy tiene 14 elementos cada uno: siempre pensé que una spline suavizada tendría tantos nudos como puntos de ajuste. Sin embargo, dado que los primeros tres y los últimos tres nudos son idénticos, 20-6 = 14, lo que tiene sentido. El problema es que no sé usarModel$fit$knot yModel$fit$coef para hacer predicciones fuera de R. Traté de echar un vistazo apredict.smooth.spline, pero sorprendentemente eso es lo que obtengo

> library(splines)
> predict.smooth.spline
Error: object 'predict.smooth.spline' not found

EDITAR: dado que aparentemente algunos usuarios malinterpretaron la pregunta, sé cómo usarpredicten R, para obtener nuevos valores de mi spline de suavizado. El problema es que quiero hacer esas predicciones en un código externo. Por lo tanto, quería echar un vistazo al código de la funciónpredict.smooth.spline, para poder intentar reproducir el algoritmo fuera de R. Por lo general, en R puede leer el código de una función simplemente ingresando su nombre (sin argumentos y sin paréntesis) en el indicador R. Pero cuando trato de hacer eso conpredict.smooth.spline, Me sale el error anterior.

EDIT2: gracias a la gran ayuda de @ r2evans, encontré la fuente parapredict método desmooth.spline. Yo (creo que) entiendo la mayor parte:

> 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)
}

Sin embargo, tengo dos dificultades:

el.Fortran() la función llama a una subrutina Fortranbvalus cuyofuente Es bastante simple. Sin embargo,bvalus a su vez llamadasbvalue que es mucho masComplicadoy llamadasinterv cuya fuente no puedo encontrar. Malas noticias:bvalue es demasiado complicado para mí entender (definitivamente no soy un experto de Fortran). Buenas noticias: el código externo que debe reproducirsepredict.smooth.spline.fit También es un código Fortran. Si lo peor llega a ser peor, podría pedirle a mi compañero de trabajo que incluya la fuente debvalus ybvalue en su código Sin embargo, incluso en este escenario ciertamente no tan agradable, todavía extrañaría el código fuente deinterv (¡Espero que no llame a otra cosa!).

No entiendo lo que se está haciendo aquí (tenga en cuenta que solo estoy interesado en elderiv == 0 caso):

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)
        }

¿Algún tipo de código recursivo? ¿Alguna ayuda aquí?