구간 \([z_i, z_{i+1}]\)에서 3차 스플라인은 3차 다항식이며 연속인 1차 및 2차 도함수를 가짐
각 매듭은 3개의 제약식을 만족
매듭이 \(n\)개인 경우 3차 스플라인을 나타내려면 \(n+4\)개의 모수가 필요
\(4(n + 1)\)개 자유도에서 \(3n\) 연속 조건을 뺌
B 스플라인 (basis splines): 임의의 스플라인 함수는 B 스플라인의 선형결함으로 표현 가능
R 함수 splines::bs()는 B 스플라인 기저 행렬을 생성하여 선형 회귀 적합에 포함될 수 있음 (회귀 스플라인)
자연 3차 스플라인 (natural cubic splines)
제한된 형태의 \(B-\)스플라인은 경계 매듭(boundary knots) 바깥에서 선형이 되도록 B 스플라인을 제약
내부 매듭 (internal knots) \(n\)개, 경계 매듭 2개: \(n+2\)개 모수
R 함수 splines::ns()는 자유도 df를 지정하여 설명변수의 경험분포의 분위수로서 내부 매듭을, 최대최소값으로 경계 매듭을 선택할 수 있음
주의: 스플라인 모형으로부터의 예측은 함수의 기저가 설명변수의 관측치에 따라 달라지므로 주의가 필요.
predict.lm()을 사용하면 새로운 기저함수 집합이 형성되고 적합된 계수가 잘못 적용됨
평활 스플라인 (smoothing splines)
\(n\)개의 순서쌍 \((x_i, y_i)\)에 대해 \[
\sum w_i\left[y_i-f\left(x_i\right)\right]^2+\lambda \int\left(f^{\prime \prime}(x)\right)^2 \mathrm{~d} x
\] 을 최소화하는 두 번 미분가능한 함수 \(f\)를 적합
\(x_i\)에 매듭이 있는 3차 스플라인
자료점을 보간하지 않음(\(\lambda > 0\))
자유도: \(\operatorname{tr}(S_{\lambda})\), \(S_{\lambda}\)는 \(\widehat{y}=S_{\lambda} y\)인 \(n \times n\) 행렬 (Green and Silverman, 1994, pp. 37-–8, Hastie, Tibshirani, and Friedman, 2009, Ch. 5)
R 함수 smooth.spline()은 \(\lambda\) 또는 자유도를 입력으로 받거나 그렇지 않으면 교차 검증을 통해 자동으로 \(\lambda\)를 선택
잔차 \[
Y - \alpha - \sum_{k \ne j} f_k (X_k)
\] 에 대해 산점도 평활자를 이용, 평활함수 \(f_j\)르 적합
모형의 선형 항들은 최소제곱법으로 적합
수렴될 때까지 반복
gam::gam()
모형 공식은 평활 스플라인과 loess를 나타내는 s(x), lo(x) 항을 허용
rock 자료
석유를 함유한 암석 12개 각각의 4개 단면에 대한 측정값(area, peri, shape—암석 단면에 있는 기공의 원형도) 포함
목표: 다른 세 가지 측정값을 통해 투과성 perm(유체 흐름의 특성)을 예측
투과율의 스케일(6.3–1300)이 크게 변하므로 로그 척도 사용
rock.lm <-lm(log(perm) ~ area + peri + shape, data = rock) rock.gam <- gam::gam(log(perm) ~s(area) +s(peri) +s(shape), control = gam::gam.control(maxit =50), data = rock)summary(rock.gam)
Call: gam::gam(formula = log(perm) ~ s(area) + s(peri) + s(shape),
data = rock, control = gam::gam.control(maxit = 50))
Deviance Residuals:
Min 1Q Median 3Q Max
-1.6848 -0.4702 0.1235 0.5430 1.2955
(Dispersion Parameter for gaussian family taken to be 0.7446)
Null Deviance: 126.9322 on 47 degrees of freedom
Residual Deviance: 26.0597 on 35.0005 degrees of freedom
AIC: 134.8981
Number of Local Scoring Iterations: NA
Anova for Parametric Effects
Df Sum Sq Mean Sq F value Pr(>F)
s(area) 1.000 14.577 14.577 19.5788 8.988e-05 ***
s(peri) 1.000 75.670 75.670 101.6321 6.856e-12 ***
s(shape) 1.000 0.413 0.413 0.5548 0.4613
Residuals 35.001 26.060 0.745
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Anova for Nonparametric Effects
Npar Df Npar F Pr(F)
(Intercept)
s(area) 3 0.34165 0.7953
s(peri) 3 0.94055 0.4315
s(shape) 3 1.43200 0.2500
anova(rock.lm, rock.gam)
Analysis of Variance Table
Model 1: log(perm) ~ area + peri + shape
Model 2: log(perm) ~ s(area) + s(peri) + s(shape)
Res.Df RSS Df Sum of Sq F Pr(>F)
1 44.000 31.949
2 35.001 26.060 8.9995 5.8891 0.8789 0.5528
rock.gam1 <- gam::gam(log(perm) ~ area + peri +s(shape), data = rock) anova(rock.lm, rock.gam1, rock.gam)
Analysis of Variance Table
Model 1: log(perm) ~ area + peri + shape
Model 2: log(perm) ~ area + peri + s(shape)
Model 3: log(perm) ~ s(area) + s(peri) + s(shape)
Res.Df RSS Df Sum of Sq F Pr(>F)
1 44.000 31.949
2 41.000 28.999 3.0000 2.9495 1.3205 0.2833
3 35.001 26.060 5.9995 2.9396 0.6581 0.6835
결과는 각 일변수함수가 하나의 선형 자유도와 세 개의 비선형 자유도를 가짐을 보여준다.
RSS를 31.95(선형 적합)에서 26.06(완전비선형적합)으로 줄이는 것은 추가 자유도 9에서는 유의하지 않음
peri에 대한 비선형 항은 유의하지 않아 보임
shape의 비선형 항만으로 RSS는 29.00
par(mfrow =c(2, 3), pty ="s") plot(rock.gam, se = T)plot(rock.gam1, se = T)
설명변수 벡터 \(\boldsymbol{X} = (x_1, \dots, X_p)\)의 차원이 높을 때, 가법 모형 \(Y=\alpha+\underset{j=1} { \stackrel{p}{\sum} } f_j (X_j)+\epsilon\)은 성분 \(X_j\)별로 상당한 자유도를 부여하기 때문에 전체 모형의 자유도가 매우 커지는 반면 변수간 교호 작용을 포함하지 않음
충분히 큰 \(M\)에 대해 \(\boldsymbol{X}\)에 대한 임의의 연속 함수를 근사할 수 있다(uniform approximation on a compact set, Diaconis and Shahshahani, 1984).
가법항 \(f_j (\boldsymbol{\alpha}_j^T \boldsymbol{X})\)은 한 방향을 제외한 모든 방향에서 일정하기 때문에 능형 함수(ridge function)로 불림
rock1 <- rock |> dplyr::mutate(area = area /10000, peri = peri /10000)(rock.ppr <-ppr(log(perm) ~ area + peri + shape, data = rock1, nterms =2, max.terms =5))
Call:
ppr(formula = log(perm) ~ area + peri + shape, data = rock1,
nterms = 2, max.terms = 5)
Goodness of fit:
2 terms 3 terms 4 terms 5 terms
8.737806 5.289517 4.745799 4.490378
summary(rock.ppr)
Call:
ppr(formula = log(perm) ~ area + peri + shape, data = rock1,
nterms = 2, max.terms = 5)
Goodness of fit:
2 terms 3 terms 4 terms 5 terms
8.737806 5.289517 4.745799 4.490378
Projection direction vectors ('alpha'):
term 1 term 2
area 0.34357179 0.37071027
peri -0.93781471 -0.61923542
shape 0.04961846 0.69218595
Coefficients of ridge terms ('beta'):
term 1 term 2
1.6079271 0.5460971
\(\beta_{ij}\): 다중 응답 회귀 모형에 대응 \[
Y_i = \alpha_{i0} + \sum_{j=1}^{M} \beta_{ij} f_j (\boldsymbol{\alpha}_j^T \boldsymbol{X}) + \epsilon
\]
능형함수 \(f_j\)는 자료의 사영에 대해 평균 0 및 분산 1이 되도록 스케일링
par(mfrow =c(3, 2))plot(rock.ppr)plot(update(rock.ppr, bass =5)) # supsmu control# use smoothing spline with GCVplot(rock.ppr2 <-update(rock.ppr, sm.method ="gcvspline", gcvpen =2))
summary(rock.ppr2)
Call:
ppr(formula = log(perm) ~ area + peri + shape, data = rock1,
nterms = 2, max.terms = 5, sm.method = "gcvspline", gcvpen = 2)
Goodness of fit:
2 terms 3 terms 4 terms 5 terms
22.52331 21.56513 21.56516 0.00000
Projection direction vectors ('alpha'):
term 1 term 2
area 0.34852171 0.44220823
peri -0.93698661 -0.85610648
shape 0.02426324 -0.26745763
Coefficients of ridge terms ('beta'):
term 1 term 2
1.4647599 0.2008702
Equivalent df for ridge terms:
term 1 term 2
2.73 2.00
2항 모형이 효과적인 것으로 보이고, area, peri가 지배적
능형함수 시각화
summary(rock1) # to find the ranges of the variables
area peri shape perm
Min. :0.1016 Min. :0.03086 Min. :0.09033 Min. : 6.30
1st Qu.:0.5305 1st Qu.:0.14149 1st Qu.:0.16226 1st Qu.: 76.45
Median :0.7487 Median :0.25362 Median :0.19886 Median : 130.50
Mean :0.7188 Mean :0.26822 Mean :0.21811 Mean : 415.45
3rd Qu.:0.8870 3rd Qu.:0.39895 3rd Qu.:0.26267 3rd Qu.: 777.50
Max. :1.2212 Max. :0.48642 Max. :0.46413 Max. :1300.00
Xp <-expand.grid(area =seq(0.1, 1.2, 0.05), peri =seq(0, 0.5, 0.02), shape =0.2)rock.grid <-cbind(Xp, fit =predict(rock.ppr2, Xp)) lattice::wireframe(fit ~ area + peri, rock.grid,screen =list(z =160, x =-60), aspect =c(1, 0.5), drape = T)
Wilcoxon signed rank test with continuity correction
data: res2^2 and res3^2
V = 3313, p-value = 0.1705
alternative hypothesis: true location shift is greater than 0
rang: Wts가 누락된 경우 runif(n,-rang, rang)의 임의의 가중치를 사용한다.
decay: 가중치붕괴 모수 \(\lambda\)
maxit: 최적화에 대한 최대 반복 횟수
Hess: 해의 헤시안 행렬이 반환되어야 하는가?
trace: 적합과정 출력?
rock 자료
rock1 <- rock |> dplyr::mutate(area = area /10000, peri = peri /10000)rock.nn <- nnet::nnet(log(perm) ~ area + peri + shape, data = rock1, size =3, decay =1e-3, linout = T, skip = T,maxit =1000, Hess = T)
# weights: 19
initial value 1519.181470
iter 10 value 32.184361
iter 20 value 21.786682
iter 30 value 16.573606
iter 40 value 14.773713
iter 50 value 14.653865
iter 60 value 14.532256
iter 70 value 14.493436
iter 80 value 14.486145
iter 90 value 14.483477
iter 100 value 14.454209
iter 110 value 14.316777
iter 120 value 13.904518
iter 130 value 13.507059
iter 140 value 13.488607
iter 150 value 13.469334
iter 160 value 13.432777
iter 170 value 13.267998
iter 180 value 13.015840
iter 190 value 12.861436
iter 200 value 12.791917
iter 210 value 12.787941
iter 220 value 12.786843
iter 230 value 12.786574
iter 240 value 12.786264
iter 250 value 12.786038
final value 12.785771
converged
