以与R的bs（）函数相同的方式在Matlab中计算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

您的代码中有两处错误

我认为 degree 和 order 之间有些混淆。您在 R 代码中正确指定了 degree=3，但在 MATLAB 中传递给 spcol 的参数是样条的 order。一般来说，n阶的样条函数是n-1阶的分段多项式。 [1]

因为 MATLAB 的 spcol 接受 order 作为输入，您需要指定 order=4 而不是您认为已经完成的 degree=3！ 您在 MATLAB 中生成了二次样条，在 R 中生成了三次样条。
看起来你的 R 图中的末端节点具有非奇异的多重性，我的意思是它们是重复的。使端点具有 degree+1 的多重性（在我们的例子中为 4）意味着它们的位置与控制多边形重合，这些被称为 clamped 结。 [2]

在 R documentation for bs 中声明结数组包含 internal 断点。看起来边界结被定义为在您的较长代码示例中被夹住，因为它们在这一行上重复 degree+1 次：
```
knots <- c(rep(Boundary.knots[1], (degree+1)), interior.knots.sorted, rep(Boundary.knots[2], (degree+1)))
```
这对于钳位端点来说是有意义的，并且支持早先关于使用 degree 输入的观点。

所以我们在 MATLAB 中的节点向量（带有夹紧端点）应该是：
```
k = [0, 0, 0, linspace(0,1,20), 1, 1, 1]
```

结论

让我们使用 order 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 样条页面

[2]：麻省理工学院的有用页面，更深入地描述了夹紧节点和数学。

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

Calculating B spline basis in Matlab in the same way as R's bs() function

r

matlab

math

spline

bspline