Почему существует одинаковое количество базисных сплайнов, заданных bs (x, intercept = FALSE) и bs (x, intercept = TRUE)? - PullRequest
Новая
       13

Почему существует одинаковое количество базисных сплайнов, заданных bs (x, intercept = FALSE) и bs (x, intercept = TRUE)?

0 голосов
/ 16 апреля 2019

Меня немного смущает функция "bs" в пакете "splines".Исходный код говорит мне, что он удалит первый столбец матрицы, когда «intercept = FALSE», но когда я попытался «intercept = FALSE», он дает другую матрицу с таким же количеством столбцов (я ожидал, что необработанная матрица, первый столбец которойне был удален).

Я сослался на «ПРАЙМЕР НА ШПИНГАХ РЕГРЕССИИ», написанный ДЖЕФФРИ С. РАСИНОМ, и сделал некоторые корректировки его кода для написания моей собственной функции «bspline».Когда я сравнил ее с функцией «bs», я обнаружил, что они возвращают ту же матрицу, когда «intercept = FALSE» (оба удаляют первый столбец необработанной матрицы), но разные матрицы, когда «intercept = TRUE». 

### my own function to define B-spline

basis <- function(x, degree, i, knots) {
  if(degree == 0){
    B <- ifelse((x >= knots[i]) & (x < knots[i+1]), 1, 0)
  } else {
    if((knots[degree+i] - knots[i]) == 0) {
      alpha1 <- 0
    } else {
      alpha1 <- (x - knots[i])/(knots[degree+i] - knots[i])
    }
    if((knots[i+degree+1] - knots[i+1]) == 0) {
      alpha2 <- 0
    } else {
      alpha2 <- (knots[i+degree+1] - x)/(knots[i+degree+1] - knots[i+1])
    }
    B <- alpha1*basis(x, (degree-1), i, knots) + alpha2*basis(x, (degree-1), (i+1), knots)
  }
  return(B)
}

bspline <- function(x, K=NULL, degree=3, interior.knots=NULL, 
                    intercept=FALSE, Boundary.knots = range(x)) {
  # The parameter K is the number of interior knots
  if(degree < 1) stop("The spline degree must be at least 1")
  if (is.null(K) && is.null(interior.knots)) stop("The parameter K and interior.knots must be given at least one")
  if (!is.null(K) && is.null(interior.knots)) {                 
    interior.knots <- if (K > 0L) {
      interior.knots <- seq.int(from = 0, to = 1, length.out = K+2L)[-c(1L, K + 2L)]  
      quantile(x, interior.knots)  
    }
  }
  if(!is.null(K) && (K<length(interior.knots))){
    warning("K<length(interior.knots)")
  }
  knots <- sort(c(rep(Boundary.knots, (degree+1)), interior.knots)) 
  n <- length(interior.knots) + degree + 1   # number of control points
  B.mat <- matrix(0,length(x),n)
  for(i in 1:n) B.mat[,i] <- basis(x, degree, i, knots)
  if(any(x == Boundary.knots[2])) B.mat[x == Boundary.knots[2], n] <- 1
  if(intercept == FALSE) {
    return(B.mat[,-1])
  } else {
    return(B.mat)
  }
}

### compare bs function and my bspline function
library("splines")
t <- sort(runif(600,0,1))
Bf <- bs(t, df=6, degree=3, intercept = FALSE)
BF <- bspline(t,K=3,degree=3,intercept = FALSE)
Bt <- bs(t, df=6, degree=3, intercept = TRUE)
BT <- bspline(t,K=3,degree=3,intercept = TRUE)

### same results when "intercept = FALSE"
Bf[1:10,1];BF[1:10,1];BT[1:10,2];

### different results when "intercept = TRUE"
Bt[1:10,1];BT[1:10,1];
dim(Bt)[2];dim(BT)[2];

Я ожидаю, что выход dim (Bt) [2] будет равен 7, но фактический вывод равен 6. На самом деле, я ожидаю, что Bt = BT.

...