在Matlab中以与R的bs()函数相同的方式计算B样条基

Question

我在Matlab中寻找(一个理想的内置)函数,以与R中相同的方式计算B样条基矩阵,例如对于具有20个等间距3度的样条的样条基础,我会在R中做

require(splines)
B = bs(x = seq(0,1,length.out=100),
        knots = seq(0, 1, length.out=20), # 20 knots
        degree = 3,
        intercept = FALSE)
matplot(B,type="l")

为了在Matlab中得到相同的结果,我想我可以使用它

B = spcol(linspace(0,1,20),3,linspace(0,1,100));
plot(B);

但是可以看出边界结然后丢失了.有任何想法,Matlab中的等效语法与R中的相同吗？

PS R用于的代码bs()有点简化:

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

bs <- function(x, degree=3, interior.knots=NULL, intercept=FALSE, Boundary.knots = c(0,1)) {
  if(missing(x)) stop("You must provide x")
  if(degree < 1) stop("The spline degree must be at least 1")
  Boundary.knots <- sort(Boundary.knots)
  interior.knots.sorted <- NULL
  if(!is.null(interior.knots)) interior.knots.sorted <- sort(interior.knots)
  knots <- c(rep(Boundary.knots[1], (degree+1)), interior.knots.sorted, rep(Boundary.knots[2], (degree+1)))
  K <- length(interior.knots) + degree + 1
  B.mat <- matrix(0,length(x),K)
  for(j in 1:K) B.mat[,j] <- basis(x, degree, j, knots)
  if(any(x == Boundary.knots[2])) B.mat[x == Boundary.knots[2], K] <- 1
  if(intercept == FALSE) {
    return(B.mat[,-1])
  } else {
    return(B.mat)
  }
}

Answer 1

您的代码中有两件事情出错了

1. 我认为程度和秩序之间存在一些混淆.您degree=3在R代码中正确指定,但在MATLAB中,传递给的参数spcol是样条曲线的顺序.通常,阶数的样条函数n是度的分段多项式n-1.[ 1 ]

   因为MATLAB spcol接受订单作为输入,所以你需要指定order=4而不是你认为你做的是什么degree=3!您在MATLAB中生成了二次样条曲线,在R中生成了三次样条曲线.

2. 看起来你的R图中的结节具有非奇异的多重性,我的意思是它们被重复.使端点具有多重性degree+1(在我们的情况下为4)意味着它们的位置与控制多边形一致,这些被称为夹紧结.[ 2 ]

   在它的R文档中bs声明knots数组包含内部断点.看起来边界节被定义为在较长的代码样本中被钳制,因为它们degree+1在这一行上重复了一次:

   knots <- c(rep(Boundary.knots[1], (degree+1)), interior.knots.sorted, rep(Boundary.knots[2], (degree+1)))
   

   这对于钳位端点是有意义的,并且支持使用度输入的早期点.

   所以在MATLAB中我们的结矢量(带夹点的终点)应该是:

   k = [0, 0, 0, linspace(0,1,20), 1, 1, 1]
   

结论

让我们使用顺序 4,并在终点处使用带结夹的结矢量:

B = spcol([0, 0, 0, linspace(0,1,20), 1, 1, 1], 4, linspace(0,1,100)); 
plot(B);

我们现在可以看到它们在R图中的边界结,以及由于3度夹紧结的影响而在每个末端的2个额外峰.

进一步阅读

[ 1 ]:维基百科B-spline页面

[ 2 ]:麻省理工学院的有用页面,更深入地描述了钳位节点和数学.