Как создать пользовательскую реализацию сплайна для передачи в формулу - PullRequest
Новая
       28

Как создать пользовательскую реализацию сплайна для передачи в формулу

1 голос
/ 03 мая 2019

Я создаю новую платформу моделирования для непараметрической регрессии.Функция pspl принимает одномерный вектор данных вместе с некоторыми дополнительными параметрами сложности и возвращает матрицу базисных функций, оцененных в данных.Цель состоит в том, чтобы эта функция использовалась точно так же, как ns или bs, т. Е. Она должна быть передана функции моделирования, такой как lm внутри формулы.Функция pspl поставляется с методом predict, который оценивает набор базисных функций для некоторых новых данных.Такая функция необходима при прогнозировании по модели, которая установлена ​​с использованием pspl.

Подгонка модели работает нормально.К сожалению, прогнозирования нового набора данных нет.Причина, по-видимому, в том, что R не находит / не использует мою функцию предиката.pspl.

Я не очень знаком с тем, как настраивать свои собственные классы и методы, поэтому некоторая помощь будет высоко ценится здесь.

library(fda)

pspl <- function(x, num.breaks = NULL, breaks = NULL, norder = 4){

  nx <- names(x)
  if(is.null(breaks)){
    breaks <- quantile(x, probs = seq(0,1,length.out = num.breaks))
    names(breaks) <- NULL
  } 
  basis <- create.bspline.basis(norder=norder, breaks=breaks)
  Z <- getbasismatrix(x, basis, nderiv=0, returnMatrix=FALSE)
  dimnames(Z) <- list(names(x), 1:ncol(Z))
  a <- list(num.breaks=num.breaks, breaks=breaks, norder=norder)
  attributes(Z) <- c(attributes(Z), a)
  class(Z) <- "pspl"
  Z
}

predict.pspl <- function(object, newx, ...){
  if(missing(newx)) 
    return(object)
  a <- c(list(x = newx), attributes(object)[c("breaks", "norder")])
  do.call("pspl", a)
}

Пример выполнения: 

set.seed(1)
x <- 2*rnorm(100)
y <- sin(x) + .5*rnorm(100)
xnew <- rnorm(100)


mod <- lm(y ~ -1 + pspl(x, num.breaks = 5))
plot(x,y)
points(x,fitted(mod), col="blue") ## correct fit
points(xnew,predict(pspl(x, num.breaks = 5), xnew)%*%coef(mod), col="green") ## what should be predicted
points(xnew, predict(mod, newdata = data.frame(x=xnew)), col="red") ## what is actually predicted

черный: данные, синий: соответствует, зеленый: правильный прогноз (сделано «вручную»), красный: фактический прогноз

...