FDA 2.0 — 탐색적 FDA의 심화 주제 개관

미분, 벌점 스무딩, 곡선 정렬: 스칼라 유사체가 없는 FDA 고유 도구

Kokoszka & Reimherr (2017) Ch.2의 세 핵심 주제를 개관한다. 함수의 미분 정보 활용, 벌점 스무딩(penalized smoothing)을 통한 노이즈 제거, 곡선 정렬(curve alignment)을 통한 위상 변동 분리를 다룬다. 각 개념이 왜 필요한지, 전통적 다변량 분석과 어떻게 다른지를 직관적으로 설명한다.

Statistics
Functional Data Analysis
저자

Kwangmin Kim

공개

2026년 04월 29일

1 이 장의 위치와 목적

Chapter 1에서는 FDA의 첫걸음을 뗐다: 기저 전개로 이산 데이터를 함수로 만들고, 표본 평균·공분산·EFPC로 요약하고, BOA 주식과 DTI 데이터에 적용했다. 이 도구들은 전통적 다변량 분석의 자연스러운 확장이다 — 벡터의 평균을 함수의 평균으로, 공분산 행렬을 공분산 함수로 바꾼 것이다.

Chapter 2는 스칼라 데이터에는 직접적 유사체가 없는 FDA 고유의 탐색적 도구를 도입한다:

도구 핵심 질문 스칼라 유사체
미분 “곡선이 얼마나 빨리 변하는가?” 없음 — 숫자에는 기울기가 없다
벌점 스무딩 “노이즈를 어떻게 제거하면서 구조를 보존하는가?” 이동평균 (매우 제한적)
곡선 정렬 “곡선들의 타이밍 차이를 어떻게 보정하는가?” 없음 — 숫자에는 위상이 없다

이 세 도구가 FDA를 단순한 “고차원 데이터 분석”과 구별짓는 핵심이다. 100차원 벡터를 다루는 것과 \([0, 1]\) 위의 함수를 다루는 것의 근본적 차이는, 함수에는 연속적 도메인에서의 순서와 매끄러움이라는 추가 구조가 있다는 점이다. 미분은 이 순서를 이용하고, 스무딩은 매끄러움을 부과하며, 정렬은 도메인의 비선형 변환을 허용한다.


2 미분 (Derivatives)

2.1 왜 미분이 가능한가

전통적 다변량 분석에서 \(\mathbf{x} = (x_1, x_2, \ldots, x_p)^T\) 의 “미분”은 정의조차 되지 않는다 — 인덱스 \(1, 2, \ldots, p\) 사이에 연속적 순서가 없기 때문이다. 하지만 함수 \(x(t)\) 는 연속 변수 \(t\) 위에 정의되므로 \(x'(t), x''(t), \ldots\) 가 자연스럽게 존재한다.

이것이 FDA의 첫 번째 고유 이점이다: 데이터의 변화율(rate of change) 정보를 직접 추출할 수 있다.

2.2 기저 전개를 이용한 미분

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

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

\(k\) 차 도함수는 기저 함수의 \(k\) 차 도함수를 이용하여 근사할 수 있다:

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

직관: 곡선의 미분을 구하려면 원시 데이터에서 수치 미분(차분)을 할 수도 있지만, 이는 노이즈를 극적으로 증폭시킨다 — 미분은 본질적으로 고주파를 강조하는 연산이기 때문이다. 기저 전개를 먼저 적용하면 고주파 노이즈가 이미 제거된 상태에서 미분하는 것이므로, 깨끗한 도함수를 얻을 수 있다.

이를 비유하면: 울퉁불퉁한 산길의 경사를 재려면, 먼저 큰 돌멩이(노이즈)를 치우고 전체 지형(신호)만 남긴 뒤 경사를 측정하는 것이 합리적이다.

2.3 기저 선택의 제약

모든 기저가 무한히 미분 가능한 것은 아니다:

기저 미분 가능 횟수 적합한 상황
Fourier 무한 주기적 데이터
B-spline (order \(p\) ) \(p - 2\) 회 연속 비주기적 데이터
Wavelet 유한 (타입에 의존) 불연속이 있는 데이터

Cubic B-spline(order 4)은 2차 도함수까지 연속이므로, 1차 미분에는 충분하다. 하지만 2차 미분(가속도)이 필요하면 order 6 이상의 B-spline을 사용해야 한다.

실무적 결론: 가속도 분석이 필요한 연구(성장 곡선, 역학 데이터 등)에서는 기저 전개 단계에서부터 높은 차수를 선택해야 한다.

2.4 미분의 응용 예: Matern 과정

Matern 과정의 매끄러움 모수 \(\nu = 5/2\) 에서의 공분산:

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

이 과정은 2회 미분 가능( \(\nu > k\) 이면 \(k\) 회 연속 미분 가능)하므로, B-spline 기반 도함수와 수치 도함수가 매우 유사하다 (Kokoszka & Reimherr, 2017, Example 2.1.1).

library(fda)

# Matern(5/2) 시뮬레이션 (100개 시점)
tgrid <- seq(0, 1, length.out = 100)

# 공분산 행렬 구성
matern52_cov <- function(h) {
  (1 + sqrt(5)*h + 5/3*h^2) * exp(-sqrt(5)*h)
}
Sigma <- outer(tgrid, tgrid, function(t, s) matern52_cov(abs(t - s)))

# 시뮬레이션
L <- chol(Sigma + 1e-10 * diag(100))
x <- as.vector(t(L) %*% rnorm(100))

# B-spline 기저로 함수 변환 (order 6 — 2차 미분을 위해)
basis <- create.bspline.basis(c(0, 1), nbasis = 50, norder = 6)
fd_obj <- smooth.basis(tgrid, x, basis)$fd

# 도함수 계산
fd_deriv <- deriv.fd(fd_obj, 1)

# 수치 도함수와 비교
numeric_deriv <- diff(x) / diff(tgrid)

# 시각화
par(mfrow = c(1, 2))
plot(fd_obj, main = "Original Curve", xlab = "t", ylab = "x(t)")
points(tgrid, x, cex = 0.5, col = "grey")

plot(fd_deriv, main = "Derivative", xlab = "t", ylab = "x'(t)")
points(tgrid[-1] - diff(tgrid)/2, numeric_deriv,
       cex = 0.5, col = "grey")
legend("topright", c("B-spline deriv", "Numeric deriv"),
       lty = c(1, NA), pch = c(NA, 1), col = c("black", "grey"))

핵심 관찰: B-spline 기반 도함수(실선)와 수치 도함수(점)가 거의 일치한다. 이는 B-spline이 원래 곡선을 충실히 표현하고 있으며, 미분 연산이 기저 함수의 미분으로 정확하게 계산됨을 보여준다.

2.5 미분이 줄 수 있는 정보

미분은 원래 곡선에서는 보이지 않는 구조를 드러낸다:

  • 1차 도함수 \(x'(t)\): 변화 속도 — “이 시점에서 곡선이 얼마나 빨리 증가/감소하는가”
  • 2차 도함수 \(x''(t)\): 가속도 — “변화의 속도가 빨라지는가 느려지는가”

성장 곡선에서 원래 곡선은 단조 증가하여 별 정보를 주지 않지만, 가속도 곡선은 사춘기 성장 급증의 시작·정점·끝을 명확히 보여준다.


3 벌점 스무딩 (Penalized Smoothing)

3.1 문제의 설정

Chapter 1의 기저 전개 \(x(t) \approx \sum_{m=1}^{M} c_m B_m(t)\) 에서 \(M\) 을 크게 잡으면 데이터의 모든 점을 지나가는 “과적합(overfitting)” 곡선이 나온다. \(M\) 을 작게 잡으면 지나치게 단순한 형태만 표현된다. 둘 다 만족스럽지 않다.

핵심 딜레마: “데이터에 충실하면서도 노이즈는 제거하는” 적절한 매끄러움을 어떻게 달성하는가?

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

관측 모형: \(y_j = x(t_j) + \varepsilon_j\) (참 곡선 + 노이즈)

벌점 스무딩은 \(K\) 개의 기저 함수(보통 관측점 수만큼 많이)를 사용하되, 곡선이 지나치게 “흔들리는(wiggly)” 것에 벌점을 부과한다:

\[ \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{매끄러움 벌점 (smoothness penalty)}} \]

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

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

3.3 직관적 이해: 적합도와 매끄러움의 줄다리기

\(\text{PSS}_\lambda\) 의 두 항은 상충한다:

  • 첫째 항 (적합도): “데이터 점에 가까이 지나가라” — 이것만 최소화하면 관측값을 정확히 통과하는 울퉁불퉁한 곡선
  • 둘째 항 (벌점): “곡선을 매끄럽게 유지하라” — 이것만 최소화하면 직선 같은 밋밋한 곡선

\(\lambda\) 는 이 줄다리기의 균형점을 결정하는 조절 모수(tuning parameter)이다:

\(\lambda\) 효과 결과
\(\lambda = 0\) 벌점 없음 데이터 점을 모두 통과하는 거친 곡선 (과적합)
\(\lambda\) 작음 약한 벌점 데이터를 비교적 충실히 따르되 약간 평활
\(\lambda\) 적절 균형 노이즈 제거, 신호 보존
\(\lambda\) 강한 벌점 거의 직선 (과평활)
\(\lambda \to \infty\) 벌점 지배 \(L(x) = 0\) 의 해 — \(L = D^2\) 이면 직선

비유하면: \(\lambda\) 는 “줄의 장력”이다. 장력 0이면 줄이 데이터 점 하나하나를 통과하며 꼬불꼬불 휘고, 장력이 무한대면 팽팽하게 직선으로 펴진다. 적절한 장력에서 줄은 데이터의 전체적 흐름을 따르면서도 개별 점의 잡음에는 반응하지 않는다.

3.4 벌점 연산자의 선택

가장 흔한 선택은 2차 미분 벌점 \(L(x)(t) = x''(t)\) 이다. 이는 “가속도가 큰 곡선을 벌한다” — 급격한 방향 전환을 억제한다.

주기적 데이터(기후, 연간 패턴)에는 조화 가속 연산자(harmonic acceleration operator) 가 적합하다:

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

이 연산자는 기본 주파수 \(\omega\) 의 sin/cos에는 벌점을 부과하지 않고, 고주파 성분에만 강한 벌점을 부과한다. 수학적으로 보면:

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

핵심 관찰: \(j = 1\) 항에는 벌점이 없다( \(j^2 - 1 = 0\) 이므로). 이는 기본 연간 주기를 있는 그대로 보존한다는 뜻이다. \(j\) 가 커질수록 \(j^2(j^2-1)^2\) 가 급격히 증가하여, 고주파 성분을 강하게 억제한다.

직관: “1년 주기의 계절 변동은 당연하니 벌하지 않지만, 하루 단위의 급변은 노이즈일 가능성이 높으니 벌한다.”

3.5 일반화 교차 검증 (GCV)

\(\lambda\) 를 어떻게 선택하는가? 이상적으로는 leave-one-out 교차 검증을 하겠지만, 매번 모형을 재적합하는 것은 계산적으로 부담이다.

벌점 스무딩은 선형 추정량(linear smoother) 이므로:

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

여기서 \(H_\lambda\)\(N \times N\) 스무딩 행렬(hat matrix)이다. 이 구조를 이용하면 교차 검증을 한 번의 행렬 연산으로 근사할 수 있다:

\[ \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} \]

직관: 분자는 잔차 제곱합(적합도), 분모는 “유효 자유도에 의한 보정”이다. \(\text{trace}(H_\lambda)\) 는 유효 모수의 수를 나타내므로, 모수가 많을수록( \(\lambda\) 가 작을수록) 분모가 작아져 GCV가 증가한다. 이것이 과적합에 대한 벌점 역할을 한다.

\(\lambda\) 를 여러 후보값에서 GCV를 계산하고 최소점을 찾으면 된다. 실무에서는 \(\log_{10} \lambda\) 를 4~9 범위에서 탐색하는 것이 일반적이다.

3.6 R 코드 예시: 캐나다 강수량 데이터

library(fda)

# 캐나다 일별 평균 로그 강수량
logprecav <- CanadianWeather$dailyAv[, , "log10precip"]

# 포화 기저 (365개 Fourier 기저 — 관측점 수와 동일)
yearRng <- c(0, 365)
daybasis <- create.fourier.basis(yearRng, nbasis = 365)

# 조화 가속 연산자 설정
Lcoef <- c(0, (2*pi/365)^2, 0, 1)
harmaccelLfd <- vec2Lfd(Lcoef, yearRng)

# GCV로 최적 lambda 탐색
loglam <- 4:9
gcvsave <- numeric(length(loglam))

for (i in seq_along(loglam)) {
  lambda <- 10^loglam[i]
  fdParobj <- fdPar(daybasis, harmaccelLfd, lambda)
  smoothlist <- smooth.basis(day.5, logprecav, fdParobj)
  gcvsave[i] <- sum(smoothlist$gcv)
}

# 최적 lambda 확인
best_idx <- which.min(gcvsave)
cat("Best log10(lambda):", loglam[best_idx], "\n")

# 최적 lambda로 평활화
lambda_opt <- 10^loglam[best_idx]
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: Raw vs Smoothed")
axisIntervals(1)
axis(2)
lines(logprec_fd[1], lwd = 3, col = "purple")  # 평활 곡선

해석: 보라색 곡선(평활)은 일별 변동을 제거하고 연간 강수 패턴만 추출한다. St. Johns는 대서양 연안 도시로, 겨울과 여름의 강수량 차이가 크지 않지만 약간의 계절 패턴을 보인다. 날카로운 일별 점프는 모두 제거되었다 — 이것이 벌점 스무딩의 효과이다.


4 곡선 정렬 (Curve Alignment / Registration)

4.1 문제: 위상 변동 (Phase Variation)

Chapter 1의 도구(평균, 공분산, EFPC)는 모두 점별(pointwise) 연산이다 — 각 시점 \(t\) 에서 곡선 값들의 평균을 구한다. 하지만 곡선들의 “사건”이 서로 다른 시점에 발생한다면, 점별 평균은 심각하게 왜곡된다.

고전적 예시: 5개의 위상 이동된 사인 곡선 \(\sin(t - \delta_n)\) 의 점별 평균. 각 곡선의 진폭은 1이지만, 평균의 진폭은 1보다 작고 지지(support)는 더 길어진다. 평균이 개별 곡선의 형태를 전혀 대변하지 못한다.

4.2 진폭 변동 vs 위상 변동

함수 표본의 변동은 두 원천으로 분해된다:

\[ \text{MSE}_{\text{total}} = \text{MSE}_{\text{amp}} + \text{MSE}_{\text{pha}} \]

변동 유형 의미 예시
진폭 변동 (amplitude) 같은 시점에서의 값 차이 “키가 큰 아이 vs 작은 아이”
위상 변동 (phase) 사건의 타이밍 차이 “일찍 성장한 아이 vs 늦게 성장한 아이”

직관: 100m 달리기를 생각하자. 진폭 변동은 “얼마나 빨리 달리는가”의 차이이고, 위상 변동은 “출발 반응 시간”의 차이이다. 기록 분석에서 출발 반응 시간을 보정하지 않으면, 실제 주행 능력에 대한 분석이 왜곡된다.

4.3 곡선 정렬의 수학적 프레임워크

관측된 곡선 \(X_n(t)\) 이 다음과 같이 분해된다고 가정한다:

\[ X_n(t) = X_n^*(h_n(t)) \]

여기서:

  • \(X_n^*(t)\) : 정렬된(등록된) 곡선 — 위상 변동이 제거된 순수 진폭 변동
  • \(h_n(t)\) : 시간 뒤틀림 함수(time warping function) — 물리적 시간 \(t\) 를 개체별 생물학적 시간으로 변환

\(h_n\) 의 조건:

  1. \(h_n(T_1) = T_1\) , \(h_n(T_2) = T_2\) (시작과 끝은 유지)
  2. \(h_n\) 은 단조 증가 (시간의 순서 보존)
  3. \(h_n\) 은 매끄러운 함수

직관: \(h_n\) 은 고무줄처럼 시간 축을 늘이거나 줄이는 변환이다. 어떤 아이의 사춘기가 12세에 시작하고 다른 아이는 14세에 시작한다면, \(h_n\) 은 두 아이의 사춘기를 같은 “생물학적 시간”으로 정렬해준다.

4.4 랜드마크 등록 (Landmark Registration)

가장 직관적인 정렬 방법이다. “중요한 사건(landmark)”이 발생하는 개별 시점 \(t_n\) 을 찾고, 이를 공통 시점 \(t_a\) (보통 표본 평균)로 매핑하는 뒤틀림 함수를 구성한다.

성장 곡선에서의 예:

  • 랜드마크: “가속도가 마지막으로 0을 지나는 시점” (= 가장 빠르게 성장하는 나이)
  • 개별 시점 \(t_n\) : 각 아이마다 다름 (10~15세)
  • 공통 시점 \(t_a\) : 전체 평균 (≈ 11.7세)
  • 뒤틀림 함수: \((0, 0)\), \((t_n, t_a)\), \((T, T)\) 을 지나는 포물선

\[ h_n(t_n) = t_a, \quad h_n(0) = 0, \quad h_n(T) = T \]

등록된 곡선:

\[ X_n^{\text{reg}}(t) = X_n(h_n^{-1}(t)) \]

이 변환 후, 모든 곡선의 가속도가 같은 시점 \(t_a\) 에서 0을 지나므로, 점별 평균이 개별 곡선의 형태를 훨씬 잘 대변한다.

4.5 연속 등록 (Continuous Registration)

랜드마크 등록의 단점은 랜드마크를 수동으로 식별해야 한다는 것이다. 연속 등록은 뒤틀림 함수 \(h_n\) 을 자동으로 찾는다 — 등록된 곡선들이 그들의 평균에 가까워지도록 \(h_n\) 을 최적화한다.

library(fda)

# 성장 곡선 준비
age <- growth$age
heightBasis <- create.bspline.basis(c(1, 18), 35, 6, age)
heightPar <- fdPar(heightBasis, 3, 10^(-0.5))
heightSmooth <- smooth.basis(age, growth$hgtf, heightPar)

# 가속도 곡선 (2차 도함수)
accelUnreg <- deriv.fd(heightSmooth$fd, 2)

# 연속 등록 (자동)
regList <- register.fd(yfd = accelUnreg)
accelReg <- regList$regfd

# 비교 시각화
par(mfrow = c(1, 2))
plot(accelUnreg, main = "Before Registration",
     xlab = "Age", ylab = "Acceleration", ylim = c(-4, 3))
lines(mean.fd(accelUnreg), lwd = 4, col = "black")

plot(accelReg, main = "After Registration",
     xlab = "Age", ylab = "Acceleration", ylim = c(-4, 3))
lines(mean.fd(accelReg), lwd = 4, col = "black")

4.6 진폭-위상 분해

등록 전후의 분산을 비교하면 위상 변동의 기여를 정량화할 수 있다:

# 뒤틀림 함수 추출
warpFunctions <- regList$warpfd

# 진폭-위상 분해
APList <- AmpPhaseDecomp(
  xfd = accelUnreg,   # 등록 전 곡선
  yfd = accelReg,     # 등록 후 곡선
  hfd = warpFunctions # 뒤틀림 함수
)

cat("MSE_amp:", APList$MS.amp, "\n")   # 진폭 변동
cat("MSE_pha:", APList$MS.pha, "\n")   # 위상 변동
cat("Phase proportion (R^2):", APList$RSQR, "\n")  # 위상 비율

성장 곡선의 결과 (Kokoszka & Reimherr, 2017):

  • \(\text{MSE}_{\text{amp}} = 5.92\)
  • \(\text{MSE}_{\text{pha}} = 3.78\)
  • 위상 비율 = 0.39 (전체 변동의 39%가 타이밍 차이)

해석: 54명 소녀의 성장 가속도 변동 중 약 39%는 순수하게 “사춘기 시기의 차이”로 설명되고, 나머지 61%만이 “성장 속도의 개인차”이다. 정렬을 적용하면 39%의 변동을 제거할 수 있어, FPCA의 효율이 크게 향상된다.

4.7 정렬의 효과: EFPC 개수 감소

정렬 전후로 90% 분산을 설명하는 데 필요한 주성분 수를 비교하면:

정렬 전 정렬 후
PC1 설명률 55% 높아짐
90% 도달 성분 수 3+ 2 이하

위상 변동이 “인위적 복잡성”을 추가하기 때문이다. 사인 곡선 5개의 예에서, 정렬 전에는 3개 PC가 필요하지만 정렬 후에는 1개로 100% 설명된다.


5 세 도구의 상호 관계

미분, 벌점 스무딩, 곡선 정렬은 독립적으로 존재하지 않고 서로 밀접하게 연결된다:

원시 데이터 (이산, 노이즈)
    ↓ 기저 전개 (Ch.1)
함수 객체
    ↓ 벌점 스무딩 → 노이즈 제거된 매끄러운 곡선
    ↓ 미분 → 변화율/가속도 곡선 추출
    ↓ 곡선 정렬 → 위상 변동 제거
    ↓
깨끗하고 정렬된 곡선 → FPCA, 회귀, 검정에 투입

실제 분석에서는 이 순서가 권장된다:

  1. 벌점 스무딩으로 노이즈를 제거한다 (성장 곡선: 반년 단위 관측 → 연속 곡선)
  2. 미분으로 관심 있는 특성을 추출한다 (성장 곡선 → 가속도 곡선)
  3. 정렬로 위상 변동을 제거한다 (사춘기 시기 보정)
  4. 정렬된 곡선에 FPCA 등 후속 분석을 적용한다

6 더 읽을거리 (Further Reading)

Chapter 2에서 소개한 세 주제는 각각 활발한 연구 분야이다:

주제 주요 참고 문헌
곡선 정렬 리뷰 Marron et al. (2015) — 최근 발전의 체계적 정리
베이지안 정렬 Earls & Hooker (2016)
미분 방정식과 FDA Ramsay et al. (2007)
함수 분류/군집화 Wang et al. (2016), Cuevas (2014)
함수 중위수·부트스트랩 Cuevas (2014)
매니폴드 위의 FDA Ettinger et al. (2016)

Chapter 2에서 다루지 않지만 중요한 FDA 연구 방향:

  • 함수 분류(classification): 관측된 곡선을 사전 정의된 그룹에 할당 (지도 학습)
  • 함수 군집화(clustering): 곡선 표본에서 자연스러운 그룹을 식별 (비지도 학습)
  • 매니폴드 위의 함수 데이터: 도메인이 뇌 표면, 혈관 등 비유클리드 공간인 경우

7 관련 주제

선행 지식

후속 주제

Subscribe

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