Hiểu về B-splines

1. B-spline là gì?

Các đường cong Spline bao gồm những đoạn đa thức kết nối với nhau một cách trơn nhẵn tại các điểm được gọi là nút. Spline được sử dụng trong các ước lượng làm trơn các đường cong hồi quy và trong các ước lượng mô hình cộng tính mở rộng (Theo TỪ ĐIỂN CÁC THUẬT NGỮ THỐNG KÊ OXFORD, Tác giả: Yadolah Dodge, Người dịch: PGS.TS Tô Cẩm Tú, Nhà xuất bản: Đại học Quốc gia Hà Nội, 2018).

Splines được đặc trưng bởi các nút chia (knots) trên 1 khoảng xác định (domain partition), bậc (order/degree).

Các đa thức Spline bậc ba được sử dụng phổ biến nhất, trong đó, các splines này phải có đạo hàm bậc nhất và bậc hai liên tục tại các nút.

B-spline (Curry-Schoenberg B-Spines) là các hàm cơ sở cho các spline có cùng thứ tự được xác định trên các nút giống nhau, có nghĩa là tất cả các hàm spline có thể có có thể được xây dựng từ sự kết hợp tuyến tính của B-splines và chỉ có một sự kết hợp độc đáo cho mỗi chức năng spline.

Lưu ý: Bậc (degree) của B-splines được định nghĩa từ 0, trong khi thứ tự (order) của B-splines được định nghĩa từ 1.

\[ order = degree +1 \]

2. Định nghĩa B-splines, đạo hàm của B-splines và KGVT \(\mathcal{S}_{k}^{\Delta \lambda}[a, b]\)

Một chút toán nhé!

Định nghĩa dưới này được tổng hợp từ De Boor, C., & De Boor, C. (1978) (page 88, 106) và Schumaker, L. (2007) (page 124).

2.1 Định nghĩa B-splines

Let the sequence of knots: ${0}=a<{1}<<{g}<b={g+1} $

The (normalized) B-spline of degree 0 (order 1), for \(i=0, \ldots, g\) \[ B_{i}^{1}(x)= \begin{cases}1 & \text { if } x \in\left[\lambda_{i}, \lambda_{i+1}\right) \\ 0 & \text { otherwise }\end{cases} \]

The (normalized) B-spline of degree 1 (order 2), for \(i=0, \ldots, g\) \[\begin{align*} B_{i}^{2}(x) & = \dfrac{x-\lambda_{i}}{\lambda_{i+1}-\lambda_{i}} B_{i}^{1}(x) +\dfrac{\lambda_{i+2}-x}{\lambda_{i+2}-\lambda_{i+1}} B_{i+1}^{1}(x)\\ & = \begin{cases} \dfrac{x-\lambda_{i}}{\lambda_{i+1}-\lambda_{i}} & \text { if } x \in\left[\lambda_{i}, \lambda_{i+1}\right) \\ \dfrac{\lambda_{i+2}-x}{\lambda_{i+2}-\lambda_{i+1}} & \text { if } x \in\left[\lambda_{i+1}, \lambda_{i+2}\right) \\ 0 & \text { otherwise }\end{cases}\\ \end{align*}\] And \[\begin{align*} B_{i+1}^{2}(x) & = \dfrac{x-\lambda_{i+1}}{\lambda_{i+2}-\lambda_{i+1}} B_{i+1}^{1}(x) +\dfrac{\lambda_{i+3}-x}{\lambda_{i+3}-\lambda_{i+2}} B_{i+2}^{1}(x)\\ & = \begin{cases} \dfrac{x-\lambda_{i+1}}{\lambda_{i+2}-\lambda_{i+1}} & \text { if } x \in \left[\lambda_{i+1}, \lambda_{i+2}\right) \\ \dfrac{\lambda_{i+3}-x}{\lambda_{i+3}-\lambda_{i+2}} & \text { if } x \in\left[\lambda_{i+2}, \lambda_{i+3}\right) \\ 0 & \text { otherwise }\end{cases}\\ \end{align*}\]

The (normalized) B-spline of degree \(k, k \in \mathbb{N},(\) order \(k+1)\) is \[ B_{i}^{k+1}(x)=\dfrac{x-\lambda_{i}}{\lambda_{i+k}-\lambda_{i}} B_{i}^{k}(x)+\dfrac{\lambda_{i+k+1}-x}{\lambda_{i+k+1}-\lambda_{i+1}} B_{i+1}^{k}(x) \]

2.2. Đạo hàm của B-splines

• The first derivatives of the B-splines

The functions \(Z_{i}^{k+1}(x)\) for \(k \geq 0, k \in \mathbb{N}\) are the first derivatives of the B-splines; \[ Z_{i}^{k+1}(x):=\dfrac{\mathrm{d}}{\mathrm{d} x} B_{i}^{k+2}(x) \] For \(k=0\) \[ Z_{i}^{1}(x)=\left\{\begin{array}{cl} \dfrac{1}{\lambda_{i+1}-\lambda_{i}} & \text { if } x \in\left[\lambda_{i}, \lambda_{i+1}\right) \\ \dfrac{-1}{\lambda_{i+2}-\lambda_{i+1}} & \text { if } x \in\left(\lambda_{i+1}, \lambda_{i+2}\right] \end{array}\right. \] For \(k \geq 1\) \[\begin{equation} Z_{i}^{k+1}(x) = (k+1)\left(\dfrac{B_{i}^{k+1}(x)}{\lambda_{i+k+1}-\lambda_{i}}-\dfrac{B_{i+1}^{k+1}(x)}{\lambda_{i+k+2}-\lambda_{i+1}}\right) = \end{equation}\]

• Functions \(Z_{i}^{k+1}(x)\) have similar properties as B-splines \(B_{i}^{k+1}(x)\).
• The second derivatives of the B-splines The functions \(S_i^{k+1}(x)\) for \(k \geq 0, k \in \mathbb{N}\) are the second derivatives of the B-splines; \[ S_i^{k+1}(x)(x):=\dfrac{\mathrm{d^2}}{\mathrm{d} x^2} B_{i}^{k+2}(x) =\dfrac{\mathrm{d}}{\mathrm{d} x} Z_{i}^{k+1}(x) \] For \(k=0\) \[ S_i^{1}(x)(x)=0 \] For \(k \geq 1\)

\[\begin{equation*} S_i^{k+1}(x)(x) = (k+1)\left(\dfrac{1}{\lambda_{i+k+1} -\lambda_{i}} \dfrac{\mathrm{d}}{\mathrm{d} x} B_{i}^{k+1}(x) -\dfrac{1}{\lambda_{i+k+2}-\lambda_{i+1}} \dfrac{\mathrm{d}}{\mathrm{d} x} B_{i+1}^{k+1}(x) \right) \end{equation*}\]

2.3.The vector space \(\mathcal{S}_{k}^{\Delta \lambda}[a, b]\)

The vector space \(\mathcal{S}_{k}^{\Delta \lambda}[a, b]\) of polynomial splines of degree \(k>0, k \in \mathbb{N}\), defined on a finite interval \(I=[a, b]\) with the sequence of knots \(\Delta \lambda=\left\{\lambda_{i}\right\}_{i=0}^{g+1}, \lambda_{0}=a<\lambda_{1}<\ldots<\lambda_{g}<b=\lambda_{g+1}\), the dimension is \[ \operatorname{dim}\left(\mathcal{S}_{k}^{\Delta \lambda}[a, b]\right)=g+k+1 \] Then adding more knots \[ \lambda_{-k}=\cdots=\lambda_{-1}=\lambda_{0}=a, \quad b=\lambda_{g+1}=\lambda_{g+2}=\cdots=\lambda_{g+k+1} \]

A basis functions \(B_{i}^{k+1}(x), i = \overline{1, g+k+1}\).\

A spline \(s_{k}(x) \in \mathcal{S}_{k}^{\Delta \lambda}[a, b]\) in \(L^{2}(I)\) has a unique representation \[ s_{k}(x)=\sum_{i=-k}^{g} b_{i} B_{i}^{k+1}(x) . \]

3. Biểu diễn B-splines: A toy example by hand!

Mục này minh họa B-plines bậc (degree) 0, 1, 2 và 3, trên khoảng \([0,6]\) và tại các nút chia 0, 1, 4, 6 để hiểu rõ về B-Splines.

Định nghĩa khoảng và nút

require(tidyverse)
lambda0 <- 0
lambda1 <- 1
lambda2 <- 4
lambda3 <- 6
lambda4 <- lambda3

B1 <- data.frame(x = seq(0, 6, by = 0.01))

3.1 B-spline degree 0 với các nút trên.

B1 <- B1 %>% mutate(B10 = ifelse(lambda0 <= x & x< lambda1, 1, 0),
                    B11 = ifelse(lambda1 <= x & x< lambda2, 1, 0),
                    B12 = ifelse(lambda2 <= x & x< lambda3, 1, 0))

plot(B1$x, B1$B10, col = "blue",cex = 1,lty = 1,
     main = "B-Spline of degree 0 ",
     ylim = c( 0, 1.2 ),
     xlab = "x", ylab = "Value")
points(B1$x, B1$B11, col = "green",cex = 1, 
       pch = 1)
points(B1$x, B1$B12, col = "pink",cex = 1, 
       pch = 20)
legend("topright",legend = c("B10", "B11", "B12"),
       lty = c(1,1,1), col = c("blue", "green", "pink"),
      cex = 1)

3.2 B-spline degree 1 với các nút trên.

#============B2=========
B2 <- B1 %>% mutate(B20 = ((x -lambda0)/(lambda1- lambda0))*B10 +  ((lambda2 - x)/(lambda2- lambda1))*B11,
                    B21 = ((x -lambda1)/(lambda2- lambda1))*B11 +  ((lambda3 - x)/(lambda3- lambda2))*B12)

plot(B2$x, B2$B20, col = "blue",cex = 1,lty = 1,
     main = "B-Spline of degree 1 ",
     ylim = c( 0, 1.2 ),
     xlab = "x", ylab = "Value")
points(B2$x, B2$B21, col = "green",cex = 1, 
       pch = 1)

legend("topright", legend = c("B20", "B21"),
       lty = c(1,1), col = c("blue", "green"),
      cex = 1)

3.3 B-spline degree 2 với các nút trên.

#============B3=========

B3 <- B2 %>% mutate(B30 = ((x -lambda0)/(lambda2- lambda1))*B20 +  
                      ((lambda3 - x)/(lambda3- lambda1))*B21)

plot(B3$x, B3$B30, col = "blue",cex = 1,lty = 1,
     main = "B-Spline of degree 2 ",
     ylim = c( 0, 1.2 ),
     xlab = "x", ylab = "Value")
legend("topright", legend = c("B30"),
       lty = c(1,1), col = c("blue"),
      cex = 1)

4. Gói splines trong R

Đây là gói lệnh cung cấp nhiều function trong ước lượng Splines, chi tiết tại https://www.rdocumentation.org/packages/splines/versions/3.6.2. Một số hàm đơn giản như dưới này:

4.1. Hàm splineDesign(): Thiết kế các B-splines trên khoảng xác định, các nút và thứ tự. Trong đó

** knots là vị trí các nút. ** x: Tập xác định của các B-splines. ** derivs: Cấp của đạo hàm của B-splines.

Ví dụ trong R:

require(splines)

## Loading required package: splines

knots <- c(1,1.8,3:5,6.5,7,8.1,9.2,10)  # 10 => 10-4 = 6 Basis splines
x <- seq(min(knots)-1, max(knots)+1, length.out = 501)
bb <- splineDesign(knots, x = x, outer.ok = TRUE)

plot(range(x), c(0,1), type = "n", xlab = "x", ylab = "",
     main =  "B-splines - sum to 1 inside inner knots")
#mtext(expression(B[j](x) *"  and "* sum(B[j](x), j == 1, 6)), adj = 0)
abline(v = knots, lty = 3, col = "light gray")
abline(v = knots[c(4,length(knots)-3)], lty = 3, col = "gray10")
lines(x, rowSums(bb), col = "gray", lwd = 2)
matlines(x, bb, ylim = c(0,1), lty = 1)

4.2. Hàm bs(): B-Spline Basis for Polynomial Splines

Thử 1 số code sẵn trong R!!

require(stats); require(graphics)
#bs(women$height, df = 5)
summary(fm1 <- lm(weight ~ bs(height, df = 5), data = women))

## 
## Call:
## lm(formula = weight ~ bs(height, df = 5), data = women)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.31764 -0.13441  0.03922  0.11096  0.35086 
## 
## Coefficients:
##                     Estimate Std. Error t value Pr(>|t|)    
## (Intercept)         114.8799     0.2167 530.146  < 2e-16 ***
## bs(height, df = 5)1   3.4657     0.4595   7.543 3.53e-05 ***
## bs(height, df = 5)2  13.0300     0.3965  32.860 1.10e-10 ***
## bs(height, df = 5)3  27.6161     0.4571  60.415 4.70e-13 ***
## bs(height, df = 5)4  40.8481     0.3866 105.669 3.09e-15 ***
## bs(height, df = 5)5  49.1296     0.3090 158.979  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2276 on 9 degrees of freedom
## Multiple R-squared:  0.9999, Adjusted R-squared:  0.9998 
## F-statistic: 1.298e+04 on 5 and 9 DF,  p-value: < 2.2e-16

## example of safe prediction
plot(women, xlab = "Height (in)", ylab = "Weight (lb)")
ht <- seq(57, 73, length.out = 200)
lines(ht, predict(fm1, data.frame(height = ht)))

## Warning in bs(height, degree = 3L, knots = c(`33.33333%` = 62.6666666666667, :
## some 'x' values beyond boundary knots may cause ill-conditioned bases

5. So sánh kết quả tự tính và hàm sẵn trong gói

Trong quá trình học, trong phần mềm R cũng thế!! Đôi khi chúng ta phải viết chi tiết các tính toán “manually” trước khi sử dụng một function để biết chính xác hàm đang tính như thế nào.

So sánh B-plines bậc (degree) 0, 1, 2 và 3, trên khoảng \([0,6]\) và tại các nút chia 0, 1, 4, 6 bằng hai phương pháp

5.1. B-Spline of degree 1 (order= 2)

knots <-c(lambda0, lambda1, lambda2, lambda3, lambda4) 
x = seq(0, 6, by = 0.01)
bb <- splineDesign(knots, x = x, outer.ok = TRUE,
                   ord = 2, derivs =0)
op <- par(mfrow = c(1, 2))
plot(B2$x, B2$B20, col = "black",cex = 1.5,lty = 1, type = "l",
     main = "By Huong ",
     ylim = c( 0, 1.2 ),
     xlab = "x", ylab = "Value")
lines(B2$x, B2$B21, col = "red",cex = 1.5)
abline(v = knots, lty = 3, col = "light gray")


plot(range(x), c(0,1), type = "n", xlab = "x", ylab = "",
     main =  "By splineDesign() function",
     ylim = c( 0, 1.2 ))
#mtext(expression(B[j](x) *"  and "* sum(B[j](x), j == 1, 6)), adj = 0)
abline(v = knots, lty = 3, col = "light gray")
abline(v = knots[c(4,length(knots)-3)], lty = 3, col = "gray10")
#lines(x, rowSums(bb), col = "gray", lwd = 2)
matlines(x, bb, ylim = c(0,1), lty = 1)

par(op)

5.1. B-Spline of degree 2 (order= 3)

#==============Compare Huong's code and package, ord = 3

knots <-c(lambda0, lambda1, lambda2, lambda3, lambda4) 
x = seq(0, 6, by = 0.01)
bb <- splineDesign(knots, x = x, outer.ok = TRUE,
                   ord = 3, derivs =0)
op <- par(mfrow = c(1, 2))
plot(B3$x, B3$B30, col = "black",cex = 1.5, lty = 1, type = "l",
     main = "By Huong ",
     ylim = c( 0, 1.2 ),
     xlab = "x", ylab = "Value")
abline(v = knots, lty = 3, col = "light gray")


plot(range(x), c(0,1), type = "n", xlab = "x", ylab = "",
     main =  "By splineDesign() function",
     ylim = c( 0, 1.2 ))
#mtext(expression(B[j](x) *"  and "* sum(B[j](x), j == 1, 6)), adj = 0)
abline(v = knots, lty = 3, col = "light gray")
abline(v = knots[c(4,length(knots)-3)], lty = 3, col = "gray10")
#lines(x, rowSums(bb), col = "gray", lwd = 2)
matlines(x, bb, ylim = c(0,1), lty = 1)

par(op)

THANK YOU!