FDA 2.1~2.2 — 미분과 벌점 스무딩

기저 전개를 이용한 도함수 계산, 벌점화 제곱합(PSS), GCV를 통한 스무딩 모수 선택

Kokoszka & Reimherr (2017) Ch.2의 첫 두 절을 상세히 다룬다. §2.1에서는 함수의 미분이 FDA에서 왜 고유한 도구인지, 기저 함수의 미분을 통해 어떻게 안정적으로 도함수를 구하는지, 기저 선택의 제약을 설명한다. §2.2에서는 벌점 스무딩의 수학적 정식화(PSS), 선형 미분 연산자의 역할, GCV를 통한 자동 모수 선택, 조화 가속 연산자의 유도를 포함한다.

Statistics
Functional Data Analysis
저자

Kwangmin Kim

공개

2026년 04월 29일

1 미분: FDA 고유의 도구

1.1 다변량 분석에는 없는 것

전통적 다변량 데이터 \(\mathbf{x} = (x_1, x_2, \ldots, x_p)^T\) 에서 인덱스 \(1, 2, \ldots, p\) 는 단순한 레이블이다. \(x_3\) 과 \(x_4\) 사이에 “연속적 변화”라는 개념이 없으므로, “\(x\) 의 미분”은 정의될 수 없다.

하지만 함수 데이터 \(x(t)\) 에서 \(t\) 는 연속 변수(시간, 파장, 위치 등)이다. \(t\) 의 미소 변화에 대한 \(x\) 의 반응 — 즉 미분 — 이 자연스럽게 정의된다. 이것이 FDA를 “고차원 다변량 분석”과 근본적으로 구별하는 첫 번째 요소이다.

1.2 미분이 주는 정보의 질적 차이

원래 곡선과 그 도함수는 서로 다른 수준의 정보를 담는다:

도함수 차수 의미 보여주는 것
\(x(t)\) 위치/수준 “현재 상태”
\(x'(t)\) 속도/변화율 “얼마나 빨리 변하는가”
\(x''(t)\) 가속도 “변화가 빨라지는가 느려지는가”

핵심 원리: 차수가 높을수록 더 미세한 구조가 드러난다. 성장 곡선 \(x(t)\) 자체는 단조 증가하여 별 정보를 주지 않지만, 가속도 \(x''(t)\) 는 사춘기 성장 급증의 시작·정점·끝을 날카롭게 보여준다. 원래 곡선에서는 보이지 않던 “사건(event)”이 도함수에서 드러나는 것이다.

비유하면: 자동차의 위치만 보면 “어딘가로 이동 중”이라는 정보뿐이다. 속도를 보면 “가속 중인지 감속 중인지” 알 수 있고, 가속도를 보면 “브레이크를 밟았는지 액셀을 밟았는지” — 즉 운전자의 의도까지 추론할 수 있다.

2 기저 전개를 이용한 미분

2.1 핵심 아이디어

곡선 \(x_n(t)\) 가 기저 전개로 표현되어 있다면:

\[ x_n(t) = \sum_{m=1}^{M} c_{nm} B_m(t) \]

\(k\) 차 도함수는 동일한 계수 \(c_{nm}\) 과 기저 함수의 \(k\) 차 도함수로 근사된다:

\[ x_n^{(k)}(t) \approx \sum_{m=1}^{M} c_{nm} B_m^{(k)}(t) \]

이것이 강력한 이유는: 원시 데이터를 다시 볼 필요 없이, 이미 구한 기저 계수만으로 도함수를 즉시 계산할 수 있다는 점이다. 또한 기저 전개 단계에서 이미 노이즈가 제거되었으므로, 미분에 의한 노이즈 증폭 문제가 자연스럽게 해결된다.

2.2 왜 수치 미분은 위험한가

원시 데이터에서 직접 수치 미분(전진/중심 차분)을 하면:

\[ x'(t_j) \approx \frac{x(t_{j+1}) - x(t_j)}{t_{j+1} - t_j} \]

이 연산은 고주파 노이즈를 극적으로 증폭한다. 수학적으로, 주파수 \(\omega\) 인 노이즈 성분 \(\varepsilon \sin(\omega t)\) 의 미분은 \(\varepsilon \omega \cos(\omega t)\) 이므로, 진폭이 \(\omega\) 배로 커진다. 고주파일수록 폭발적으로 증폭된다.

기저 전개 후 미분은 이 문제를 우회한다: \(M\) 개의 기저 함수가 이미 대역 제한(band-limiting)을 수행하여 고주파 노이즈가 제거된 상태에서 미분하기 때문이다.

방법 노이즈 처리 결과 품질
수치 차분 증폭 거칠고 불안정
기저 전개 후 미분 이미 제거됨 매끄럽고 안정적

2.3 기저 선택의 제약 조건

모든 기저가 임의 차수의 미분을 허용하는 것은 아니다:

미분 가능 횟수와 기저 선택

기저 함수 \(B_m(t)\) 이 \(k\) 차 도함수를 가지려면, \(B_m\) 자체가 \(k\) 회 이상 연속 미분 가능해야 한다.

기저 유형 차수(order) 연속 미분 가능 횟수 1차 미분 2차 미분
Fourier - \(\infty\) 가능 가능
B-spline (cubic, order 4) 4 2 (= order - 2) 가능 경계 주의
B-spline (order 6) 6 4 가능 가능
B-spline (order 8) 8 6 가능 가능

실무적 규칙: \(k\) 차 도함수가 필요하면 order \(\geq k + 2\) 인 B-spline을 사용한다. 2차 도함수(가속도)가 목적이면 최소 order 6을 선택해야 한다.

Fourier 기저의 미분은 특히 깔끔하다:

\[ \frac{d}{dt} \sin(\omega j t) = \omega j \cos(\omega j t), \quad \frac{d}{dt} \cos(\omega j t) = -\omega j \sin(\omega j t) \]

미분 후에도 같은 종류의 삼각함수이므로, Fourier 기저는 미분에 대해 닫혀 있다(closed under differentiation). B-spline의 도함수는 한 차수 낮은 B-spline이 된다 — cubic B-spline의 도함수는 quadratic B-spline이다.

3 예시: Matern 과정의 미분 (Example 2.1.1)

3.1 Matern 과정이란

Matern 과정은 매끄러움 모수 \(\nu\) 로 곡선의 미분 가능 횟수를 제어할 수 있는 가우시안 과정이다. \(\nu > k\) 이면 \(k\) 회 연속 미분 가능한 경로를 생성한다.

\(\nu = 5/2\) , \(\sigma^2 = 1\) , \(\rho = 1\) 일 때 공분산 함수:

\[ c(t, s) = \left(1 + \sqrt{5}|t-s| + \frac{5}{3}|t-s|^2 \right) \exp(-\sqrt{5}|t-s|) \]

왜 Matern을 선택하는가? 브라운 운동( \(\nu = 1/2\) )은 미분 불가능한 거친 경로를 생성하므로 미분 예시로 부적합하다. \(\nu = 5/2\) 는 2회 미분 가능한 매끄러운 경로를 보장하여, B-spline 기반 도함수의 정확성을 검증할 수 있다.

3.2 R 코드: 시뮬레이션과 도함수 비교

library(fda)

# --- Matern(5/2) 공분산 행렬 구성 ---
matern52 <- function(h) {
  # nu = 5/2, sigma^2 = 1, rho = 1
  (1 + sqrt(5)*h + 5/3*h^2) * exp(-sqrt(5)*h)
}

n_t <- 100
tgrid <- seq(0, 1, length.out = n_t)

Sigma <- outer(tgrid, tgrid, function(t, s) matern52(abs(t - s)))

# --- 시뮬레이션: 1개 경로 생성 ---
set.seed(42)
L <- chol(Sigma + 1e-10 * diag(n_t))
x_obs <- as.vector(t(L) %*% rnorm(n_t))

# --- B-spline 기저로 함수 변환 ---
# order 6 사용: 2차 도함수까지 안전하게 계산하기 위함
basis <- create.bspline.basis(c(0, 1), nbasis = 50, norder = 6)
fd_obj <- smooth.basis(tgrid, x_obs, basis)$fd

# --- 도함수 계산: B-spline 기반 ---
fd_deriv <- deriv.fd(fd_obj, 1)  # 1차 도함수

# --- 수치 도함수: 중심 차분 ---
dt <- diff(tgrid)
numeric_deriv <- diff(x_obs) / dt
t_mid <- tgrid[-n_t] + dt/2  # 중간점

# --- 시각화 ---
par(mfrow = c(1, 2))

# 왼쪽: 원래 곡선
plot(fd_obj, main = "Matern(5/2) Realization",
     xlab = "t", ylab = "x(t)", lwd = 2, col = "red")
points(tgrid, x_obs, cex = 0.5, col = "grey40")
legend("topright", c("B-spline fit", "Observations"),
       lty = c(1, NA), pch = c(NA, 1), col = c("red", "grey40"))

# 오른쪽: 도함수 비교
plot(fd_deriv, main = "Derivative Comparison",
     xlab = "t", ylab = "x'(t)", lwd = 2, col = "red")
points(t_mid, numeric_deriv, cex = 0.5, col = "grey40")
legend("topright", c("B-spline deriv", "Numeric deriv"),
       lty = c(1, NA), pch = c(NA, 1), col = c("red", "grey40"))

3.3 결과 해석

B-spline 기반 도함수(빨간 실선)와 수치 도함수(회색 점)가 거의 완벽하게 일치한다. 이는 두 가지를 확인해준다:

  1. B-spline 적합의 충실성: 50개 기저가 원래 곡선을 거의 완벽하게 재현한다
  2. 해석적 미분의 정확성: 기저 함수의 도함수를 이용한 계산이 수치적 방법과 동치이다

Matern(5/2) 과정은 2회 미분 가능하므로 이 일치가 보장된다. 만약 \(\nu = 1/2\) (브라운 운동)에서 같은 실험을 하면, 경로가 미분 불가능하여 수치 도함수도 B-spline 도함수도 “진짜 도함수”의 근사가 아닌 인공물(artifact)이 된다.

4 벌점 스무딩의 동기

4.1 기저 전개만으로 충분하지 않은 경우

Chapter 1의 기저 전개 \(x(t) \approx \sum_{m=1}^{M} c_m B_m(t)\) 는 원시 데이터가 이미 비교적 매끄러운 경우에 잘 작동한다. 기저 수 \(M\) 을 적절히 선택하면 곡선의 전체 형태를 포착하면서 약간의 평활화가 자동으로 이루어진다.

하지만 원시 데이터에 상당한 노이즈가 있으면 문제가 생긴다:

  • \(M\) 이 작으면: 가능한 형태가 \(M\) 개 기저의 선형 결합으로 제한되어, 복잡한 구조를 놓친다
  • \(M\) 이 크면: 기저 전개가 데이터 점을 거의 정확히 통과하므로 노이즈까지 함께 “학습”한다

딜레마: 많은 기저로 풍부한 표현력을 확보하면서도, 노이즈는 제거하고 싶다. 벌점 스무딩은 이 딜레마를 정확히 해결한다.

4.2 관측 모형

단일 곡선 \(x(t)\) 를 시점 \(t_j\) 에서 노이즈와 함께 관측한다:

\[ y_j = x(t_j) + \varepsilon_j, \quad E[\varepsilon_j] = 0 \]

여기서 \(x(t)\) 는 매끄러운 참 곡선, \(\varepsilon_j\) 는 측정 오차이다. 목표는 \(\{(t_j, y_j)\}\) 로부터 \(x(t)\) 를 추정하는 것이다.

비유: 안개 낀 날 산 능선의 사진을 찍었다. 사진(데이터 \(y_j\) )에는 안개(노이즈 \(\varepsilon_j\) )와 능선(참 곡선 \(x(t)\) )이 섞여 있다. 벌점 스무딩은 “산 능선은 매끄러워야 한다”는 사전 지식(벌점)을 이용하여 안개를 걷어내는 것이다.

5 벌점화 제곱합 (Penalized Sum of Squares)

5.1 수학적 정식화

\(K\) 개의 기저 함수(종종 관측점 수와 같거나 더 많이)를 사용하여 \(x(t) \approx x_K(t) = \sum_{k=1}^{K} c_k B_k(t)\) 로 근사한다. 벌점화 제곱합(PSS)을 최소화하는 계수 \(c_1, \ldots, c_K\) 를 구한다:

\[ \text{PSS}_\lambda(c_1, \ldots, c_K) = \underbrace{\sum_j (y_j - x_K(t_j))^2}_{\text{적합도 항 (fidelity)}} + \underbrace{\lambda \int_{T_1}^{T_2} [L(x_K)(t)]^2 \, dt}_{\text{거칠기 벌점 (roughness penalty)}} \]

여기서 \(L\) 은 선형 미분 연산자:

\[ L(x)(t) = \alpha_0(t) x(t) + \alpha_1(t) x^{(1)}(t) + \cdots + \alpha_m(t) x^{(m)}(t) \]

5.2 각 항의 역할

적합도 항 \(\sum_j (y_j - x_K(t_j))^2\):

  • “추정 곡선이 데이터에 가까이 지나가도록 하라”
  • 이것만 최소화하면 \(K\) 가 충분히 크면 \(x_K(t_j) = y_j\) (데이터를 정확히 통과 = 과적합)

거칠기 벌점 \(\lambda \int [L(x_K)(t)]^2 dt\):

  • “곡선이 너무 흔들리면 벌을 주라”
  • \(L = D^2\) (2차 미분)이면: “가속도가 큰 곡선은 거칠다” → 급격한 방향 전환 억제
  • 이것만 최소화하면 \(L(x) = 0\) 의 해: \(D^2 x = 0\) 이면 \(x(t) = at + b\) (직선)

조절 모수 \(\lambda\): 두 항의 상대적 중요도를 결정한다.

5.3 \(\lambda\) 의 효과: 과적합-과평활 스펙트럼

\(\lambda\) 를 연속적으로 변화시키면 추정 곡선이 다음과 같이 변한다:

\[ \lambda = 0 \quad \xrightarrow{\text{거칠어짐}} \quad \lambda \text{ 작음} \quad \xrightarrow{\text{균형점}} \quad \lambda \text{ 적절} \quad \xrightarrow{\text{밋밋해짐}} \quad \lambda \to \infty \]

구체적으로:

\(\lambda\) 추정 곡선의 특성 극한 형태
\(0\) 모든 데이터 점을 통과 — 과적합 보간(interpolation)
\(10^2\) 약간의 평활, 대부분의 구조 보존 -
\(10^6\) (적절) 노이즈 제거, 신호의 주요 구조 보존 -
\(10^{12}\) 거의 모든 구조가 사라짐 \(L(x) = 0\) 의 해
\(\to \infty\) 벌점이 지배 \(D^2\) 이면 직선, 조화가속이면 기본 주기

직관: \(\lambda\) 는 “곡선에 걸린 고무줄의 장력”이다. 장력이 0이면 곡선이 데이터 점마다 꼬불꼬불 휘고, 장력이 무한대면 팽팽하게 직선(또는 기본 주기)으로 펴진다.

6 선형 미분 연산자 \(L\) 의 설계

6.1 가장 흔한 선택: 2차 미분 벌점

\[ L(x)(t) = x''(t) = D^2 x(t) \]

이 선택의 의미: “곡률(curvature)이 큰 곡선을 벌한다.” 결과적으로 추정 곡선은 가능한 한 곡률이 작은 — 즉 천천히 방향을 바꾸는 — 형태가 된다.

\(\lambda \to \infty\) 일 때 극한 형태는 \(x''(t) = 0\) 의 해, 즉 직선 \(x(t) = at + b\) 이다.

적합한 상황: 비주기적이고 전반적으로 매끄러운 곡선 (예: 성장 곡선의 원시 관측)

6.2 주기적 데이터를 위한 설계: 조화 가속 연산자

기후 데이터(일별 기온, 강수량)는 연간 주기 \(T\) 를 가진다. 이런 데이터에 단순 \(D^2\) 벌점을 적용하면 연간 주기 자체를 “거친 것”으로 판단하여 억제할 수 있다 — 이는 바람직하지 않다.

조화 가속 연산자(harmonic acceleration operator) 는 기본 주기를 벌하지 않도록 설계된다:

\[ L(x)(t) = \omega^2 x^{(1)}(t) + x^{(3)}(t), \quad \omega = \frac{2\pi}{T} \]

6.3 조화 가속 연산자의 핵심 성질: 유도

Fourier 전개를 사용한다:

\[ x_J(t) = c_0 + \sum_{j=1}^{J} [a_j \sin(\omega j t) + b_j \cos(\omega j t)] \]

\(L\) 을 \(\sin(\omega j t)\) 에 적용:

\[ \begin{aligned} x^{(1)}(t) &= \omega j \cos(\omega j t) \\ x^{(3)}(t) &= -\omega^3 j^3 \cos(\omega j t) \\ L(x)(t) &= \omega^2 \cdot \omega j \cos(\omega j t) + (-\omega^3 j^3 \cos(\omega j t)) \\ &= \omega^3 j (1 - j^2) \cos(\omega j t) \end{aligned} \]

핵심 관찰: \(j = 1\) 이면 \(L(\sin(\omega t)) = \omega^3 \cdot 1 \cdot (1 - 1) \cdot \cos(\omega t) = 0\)

즉, 기본 주파수 \(\omega\) 의 성분에는 벌점이 0이다. 마찬가지로 \(\cos(\omega t)\) 에도 벌점이 0이다. 상수항 \(c_0\) 에도 벌점이 없다.

6.4 벌점 적분의 계산

\(L(x_J)(t)\) 의 제곱을 \([0, T]\) 에서 적분하면:

\[ \int_0^T [L(x_J)(t)]^2 \, dt = \pi \omega^5 \sum_{j=2}^{J} j^2(j^2 - 1)^2 (a_j^2 + b_j^2) \]

이 식의 구조를 분석하면:

주파수 배수 \(j\) 벌점 가중치 \(j^2(j^2-1)^2\) 해석
1 \(1 \cdot (1-1)^2 = 0\) 기본 주기 — 벌점 없음
2 \(4 \cdot 9 = 36\) 2배 주파수 — 약한 벌점
3 \(9 \cdot 64 = 576\) 3배 주파수 — 중간 벌점
5 \(25 \cdot 576 = 14400\) 5배 주파수 — 강한 벌점
10 \(100 \cdot 9801 = 980100\) 10배 주파수 — 매우 강한 벌점

직관: 이 연산자는 “1년 주기의 계절 변동은 있는 그대로 두되, 주간·일간 변동은 점점 강하게 억제하는” 선택적 저주파 통과 필터(selective low-pass filter)이다. \(j\) 가 커질수록 벌점 가중치가 \(O(j^6)\) 으로 폭증하여, 고주파 성분의 계수 \(a_j, b_j\) 를 강제로 0에 가깝게 만든다.

6.5 \(D^2\) 벌점 vs 조화 가속 벌점 비교

\(L = D^2\) \(L = \omega^2 D + D^3\)
벌점 없는 형태 직선 \(at + b\) \(c_0 + a_1\sin\omega t + b_1\cos\omega t\)
\(\lambda \to \infty\) 극한 직선으로 수렴 기본 주기 곡선으로 수렴
적합 대상 비주기 곡선 주기 \(T\) 의 곡선
고주파 억제 강도 \(\omega_j^4\) \(j^2(j^2-1)^2 \sim j^6\)

7 일반화 교차 검증 (GCV)

7.1 모수 선택의 원칙

\(\lambda\) 가 너무 작으면 과적합, 너무 크면 과평활이다. 이상적 \(\lambda\) 는 예측 오차를 최소화하는 값이다. Leave-one-out 교차 검증(LOO-CV)은 원칙적으로 이를 달성하지만, 관측점이 \(N\) 개이면 모형을 \(N\) 번 재적합해야 하므로 계산 부담이 크다.

7.2 선형 추정량의 구조

벌점 스무딩은 선형 추정량(linear smoother) 이다: 관측 벡터 \(\mathbf{Y} = (y_1, \ldots, y_N)^T\) 에 대해 예측값이 행렬 곱으로 표현된다:

\[ \hat{\mathbf{Y}} = H_\lambda \mathbf{Y} \]

여기서 \(H_\lambda\) 는 \(N \times N\) 스무딩 행렬(hat matrix) 이다.

왜 선형인가? PSS를 \(c_k\) 에 대해 미분하고 0으로 놓으면, 적합도 항은 \(c_k\) 에 대해 이차이고 벌점 항도 이차이므로, 최적 계수 \(\hat{c}\) 는 \(\mathbf{Y}\) 의 선형 함수가 된다. 따라서 예측값 \(\hat{\mathbf{Y}}\) 도 \(\mathbf{Y}\) 의 선형 함수이다.

7.3 GCV 공식

LOO-CV를 모형 재적합 없이 근사하는 것이 GCV의 핵심이다:

\[ \text{GCV}(\lambda) = \frac{N^{-1} \sum_{n=1}^{N} (Y_n - \hat{Y}_n)^2}{\left(1 - N^{-1} \text{trace}(H_\lambda) \right)^2} \]

7.4 각 구성 요소의 해석

분자: \(N^{-1} \|\mathbf{Y} - \hat{\mathbf{Y}}\|^2\) = 평균 잔차 제곱 (적합도)

  • \(\lambda\) 가 작을수록 분자가 작아진다 (과적합이면 잔차 ≈ 0)
  • \(\lambda\) 가 클수록 분자가 커진다 (과평활이면 데이터를 못 따라감)

분모: \((1 - N^{-1} \text{trace}(H_\lambda))^2\)

  • \(\text{trace}(H_\lambda)\) = 유효 자유도(effective degrees of freedom) = 유효 모수 수
  • \(\lambda = 0\) 이면 \(H_\lambda = I\) (보간), \(\text{trace} = N\), 분모 → 0, GCV → \(\infty\)
  • \(\lambda \to \infty\) 이면 \(\text{trace} \to\) 작은 수 (2 또는 3), 분모 → 1에 가까움

분모의 역할: 과적합에 대한 자동 벌점이다. 모수를 많이 쓸수록( \(\text{trace}\) 가 클수록) 분모가 작아져 GCV가 증가한다. 이는 잔차가 작더라도 “그건 단지 데이터를 외운 것이다”라고 경고하는 메커니즘이다.

7.5 GCV vs LOO-CV

GCV가 LOO-CV를 근사하는 수학적 근거: 선형 추정량에서 \(i\) 번째 관측을 제외한 예측값은:

\[ \hat{Y}_{(-i)} = \hat{Y}_i - \frac{H_{ii}}{1 - H_{ii}} (Y_i - \hat{Y}_i) \]

LOO-CV = \(N^{-1} \sum (Y_i - \hat{Y}_{(-i)})^2 = N^{-1} \sum \frac{(Y_i - \hat{Y}_i)^2}{(1 - H_{ii})^2}\)

GCV는 개별 \(H_{ii}\) 를 평균 \(N^{-1}\text{trace}(H)\) 로 대체한 근사이다. 이 대체는 계산을 극적으로 간소화한다 — \(N\) 번의 모형 재적합이 한 번의 행렬 trace 계산으로 줄어든다.

8 실습: 캐나다 일별 강수량 데이터 (Example 2.2.2)

8.1 데이터 설명

캐나다 35개 기상 관측소의 일별 평균 로그 강수량 데이터이다. fda 패키지의 CanadianWeather 객체에 포함되어 있다. 365일 × 35개 관측소의 행렬이며, 각 열이 하나의 곡선이다.

이 데이터에 벌점 스무딩이 필요한 이유: 수십 년간 평균하더라도 일별 변동이 상당히 남아 있다. 예를 들어 6월 30일과 7월 1일의 평균 강수량이 눈에 띄게 다를 수 있지만, 이는 기후적 실체가 아닌 표본 변동성(sample variability)이다. 기후학적으로 의미 있는 것은 연간 패턴이다.

8.2 Step 1: 포화 기저로 함수 변환 (스무딩 없이)

library(fda)

# 데이터 준비
logprecav <- CanadianWeather$dailyAv[, , "log10precip"]

# 포화 기저: 관측점 수(365)와 동일한 Fourier 기저
# → 데이터를 정확히 재현 (보간), 평활화 효과 없음
yearRng <- c(0, 365)
daybasis <- create.fourier.basis(yearRng, nbasis = 365)

# 포화 기저로 변환 (스무딩 없음)
dayprecfd <- smooth.basis(
  argvals = day.5,           # 0.5, 1.5, ..., 364.5
  y = logprecav,
  fdParobj = daybasis
)$fd

# St. Johns (첫 번째 관측소) 시각화
plot(logprecav[, 1], axes = FALSE,
     xlab = "Day", ylab = "log10(mm)",
     main = "St. Johns: No Smoothing (Saturated Basis)")
lines(dayprecfd[1], col = "red", lwd = 2)
axisIntervals(1)
axis(2)

관찰: 빨간 선(포화 기저 적합)이 모든 데이터 점을 통과한다. 일별 변동이 그대로 남아 있어 연간 패턴을 파악하기 어렵다. 이것이 “과적합” 상태이다.

8.3 Step 2: 조화 가속 연산자 설정

# 조화 가속 연산자: L(x) = omega^2 * x' + x'''
# omega = 2*pi / 365
Lcoef <- c(0, (2*pi/diff(yearRng))^2, 0, 1)
# Lcoef = c(alpha_0, alpha_1, alpha_2, alpha_3)
# = c(0, omega^2, 0, 1)
# → L(x) = 0*x + omega^2*x' + 0*x'' + 1*x'''

harmaccelLfd <- vec2Lfd(Lcoef, yearRng)

코드 해석: vec2Lfd는 계수 벡터 \([\alpha_0, \alpha_1, \alpha_2, \alpha_3]\) 로부터 \(L(x) = \alpha_0 x + \alpha_1 x' + \alpha_2 x'' + \alpha_3 x'''\) 을 구성한다. 여기서 \(\alpha_1 = \omega^2, \alpha_3 = 1\) 이므로 \(L(x) = \omega^2 x' + x'''\) — 바로 조화 가속 연산자이다.

8.4 Step 3: GCV로 최적 \(\lambda\) 탐색

# lambda 후보: 10^4 ~ 10^9
loglam <- 4:9
nlam <- length(loglam)
gcvsave <- numeric(nlam)
dfsave <- numeric(nlam)

for (ilam in 1:nlam) {
  lambda <- 10^loglam[ilam]
  
  # fdPar: 기저 + 벌점 연산자 + lambda를 하나의 객체로 포장
  fdParobj <- fdPar(daybasis, harmaccelLfd, lambda)
  
  # 벌점 스무딩 적합
  smoothlist <- smooth.basis(day.5, logprecav, fdParobj)
  
  # GCV 점수와 유효 자유도 저장
  gcvsave[ilam] <- sum(smoothlist$gcv)
  dfsave[ilam] <- smoothlist$df
}

# GCV 최소점 확인
plot(loglam, gcvsave, type = "b", lwd = 2,
     xlab = expression(log[10](lambda)),
     ylab = "GCV Score",
     main = "GCV vs Smoothing Parameter")
abline(v = loglam[which.min(gcvsave)], lty = 2, col = "red")

cat("Optimal log10(lambda):", loglam[which.min(gcvsave)], "\n")
cat("Effective df at optimum:", dfsave[which.min(gcvsave)], "\n")

예상 결과: \(\log_{10} \lambda = 6\) 에서 GCV가 최소화된다 (Kokoszka & Reimherr, 2017). 이는 \(\lambda = 10^6\) 이라는 큰 값으로, 벌점이 상당히 강하게 작용하여 매끄러운 연간 패턴만 남긴다는 뜻이다.

8.5 Step 4: 최적 \(\lambda\) 로 평활 곡선 생성

# 최적 lambda 적용
lambda_opt <- 1e6
fdParobj <- fdPar(daybasis, harmaccelLfd, lambda_opt)
logprec_fd <- smooth.basis(day.5, logprecav, fdParobj)$fd

# St. Johns: 원시 데이터 + 포화 기저 + 평활 곡선 비교
plot(logprecav[, 1], axes = FALSE,
     xlab = "Day", ylab = "log10(mm)",
     main = "St. Johns: Smoothing Comparison", col = "grey60")
lines(dayprecfd[1], col = "red", lwd = 1)      # 포화 기저 (과적합)
lines(logprec_fd[1], col = "purple", lwd = 3)   # 벌점 스무딩
axisIntervals(1)
axis(2)
legend("topright",
       c("Raw data", "No smoothing", "Penalized (GCV)"),
       pch = c(1, NA, NA), lty = c(NA, 1, 1),
       col = c("grey60", "red", "purple"), lwd = c(NA, 1, 3))

8.6 결과 해석

  • 회색 점 (원시 데이터): 일별 변동이 심하여 연간 패턴을 파악하기 어렵다
  • 빨간 선 (포화 기저): 모든 점을 통과하여 과적합 — 날카로운 일별 점프가 그대로 보인다
  • 보라 선 (벌점 스무딩): 매끄러운 연간 강수 패턴만 추출 — St. Johns의 기후학적 특성을 보여준다

St. Johns는 캐나다 대서양 연안(뉴펀들랜드)에 위치하여, 전형적인 해양성 기후를 보인다. 평활 곡선은 겨울과 여름의 강수량 차이가 크지 않음을 보여준다 — 이것이 기후학적으로 의미 있는 신호이다.

8.7 다수 관측소 비교

# 첫 5개 관측소의 평활 곡선 비교
par(mfrow = c(2, 3))
for (i in 1:5) {
  plot(logprecav[, i], axes = FALSE, col = "grey60", cex = 0.3,
       xlab = "Day", ylab = "log10(mm)",
       main = CanadianWeather$place[i])
  lines(logprec_fd[i], lwd = 3, col = "purple")
  axisIntervals(1)
  axis(2)
}

# 전체 35개 평활 곡선 겹쳐 그리기
plot(logprec_fd, main = "35 Canadian Stations (Smoothed)",
     xlab = "Day", ylab = "log10(mm)", col = "grey50")
lines(mean.fd(logprec_fd), lwd = 3, col = "black")

관찰할 수 있는 패턴: - 태평양 연안 (밴쿠버): 겨울에 강수량 높음, 여름에 낮음 - 대서양 연안 (St. Johns): 연중 비교적 균일 - 내륙 (에드먼턴): 여름에 높고 겨울에 낮음

이런 기후학적 구분이 벌점 스무딩 후에야 명확히 드러난다 — 원시 데이터의 일별 노이즈가 이 구조를 가리고 있었기 때문이다.

9 유효 자유도의 해석

\(\text{trace}(H_\lambda)\) 는 “모형이 실질적으로 사용하는 모수의 수”를 나타내며, 유효 자유도(effective degrees of freedom, df) 라 부른다:

\(\lambda\) \(\text{trace}(H_\lambda)\) 해석
\(0\) \(K\) (= 기저 수) 모든 기저 사용 = 보간
\(10^6\) (최적) ~12 (예시) 12개 모수 분량의 복잡도
\(\to \infty\) 2~3 직선/기본 주기만

유효 자유도가 12라면, 365개 Fourier 기저를 사용하더라도 실질적으로는 12개 모수만 “활성화”된 셈이다. 나머지 353개 기저의 기여는 벌점에 의해 사실상 0으로 억제된다. 이것이 벌점 스무딩의 자동 모형 선택 효과이다 — 적절한 복잡도를 데이터 기반으로 결정한다.

10 PSS 최적화의 해석적 해

PSS를 최소화하는 계수 \(\hat{\mathbf{c}} = (\hat{c}_1, \ldots, \hat{c}_K)^T\) 를 해석적으로 구할 수 있다.

\(\Phi\) 를 기저 행렬 \([\Phi]_{jk} = B_k(t_j)\) , \(R\) 을 벌점 행렬 \([R]_{kk'} = \int B_k^{(m)}(t) B_{k'}^{(m)}(t) \, dt\) 라 하면:

\[ \hat{\mathbf{c}} = (\Phi^T \Phi + \lambda R)^{-1} \Phi^T \mathbf{Y} \]

이 식의 구조:

  • \(\Phi^T \Phi\) : 적합도에서 오는 항 (데이터 충실)
  • \(\lambda R\) : 벌점에서 오는 항 (매끄러움 강제)
  • \(\lambda\) 가 크면 \(R\) 의 기여가 커져 고주파 계수가 억제됨

스무딩 행렬은:

\[ H_\lambda = \Phi (\Phi^T \Phi + \lambda R)^{-1} \Phi^T \]

이므로 \(\hat{\mathbf{Y}} = H_\lambda \mathbf{Y}\) 가 성립하며, GCV 공식에 직접 대입할 수 있다.

11 벌점 스무딩의 실무적 권장사항

실무 가이드라인
  1. 기저 수를 넉넉히 잡아라: \(K \geq N\) (관측점 수 이상)도 괜찮다. 벌점이 과적합을 방지한다
  2. GCV 결과를 맹신하지 말라: 자동 선택된 \(\lambda\) 로 적합한 곡선을 반드시 시각적으로 검토하라
  3. 연산자 \(L\) 을 데이터에 맞게 설계하라: 주기 데이터에 \(D^2\) 를 쓰면 기본 주기가 과평활된다
  4. 미분과 결합하라: 벌점 스무딩 → 도함수 추출의 순서가 노이즈 제거와 미분을 동시에 달성한다

12 정리: §2.1과 §2.2의 연결

미분과 벌점 스무딩은 독립적 주제가 아니라 밀접하게 연결된다:

원시 데이터 (노이즈 포함)
    ↓
기저 전개 (기본적 차원 축소)
    ↓
벌점 스무딩 (노이즈 제거 — L이 미분 연산자)
    ↓
도함수 추출 (깨끗한 곡선에서 안정적 미분)
    ↓
후속 분석 (FPCA, 회귀, 정렬 등)

벌점 연산자 \(L\) 자체가 미분으로 정의되어 있다는 점에 주목하라. \(L = D^2\) 이면 “2차 도함수(가속도)가 큰 곡선을 벌하라”이고, 조화 가속이면 “3차 도함수와 1차 도함수의 특정 결합이 큰 곡선을 벌하라”이다. 즉, 스무딩은 미분 구조를 이용하여 매끄러움을 정의하고 부과하는 것이다.

이것이 FDA를 커널 스무딩이나 이동 평균 같은 비모수적 방법과 구별한다: FDA의 벌점 스무딩은 미분 연산자를 통해 곡선의 물리적/수학적 구조를 반영한 스무딩을 수행한다.

13 관련 주제

선행 지식

후속 주제

Subscribe

Enjoy this blog? Get notified of new posts by email: