Chapter 2 연습문제 풀이

1 문제 2.1: 등식 (2.2) 증명

1.1 문제

조화 가속 연산자(harmonic acceleration operator)에 대해 다음 등식을 증명하라:

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

여기서 \(L(x)(t) = \omega^2 x^{(1)}(t) + x^{(3)}(t)\), \(\omega = 2\pi / T\), 그리고 \(x_J(t) = c_0 + \sum_{j=1}^J [a_j \sin(\omega j t) + b_j \cos(\omega j t)]\) 이다.

1.2 풀이

Step 1: 미분 계산

푸리에 전개의 각 항에 대해 미분을 계산한다:

\[ \begin{aligned} x^{(1)}(t) &= \sum_{j=1}^J \omega j \bigl[a_j \cos(\omega j t) - b_j \sin(\omega j t)\bigr] \\ x^{(3)}(t) &= \sum_{j=1}^J (-\omega j)^3 \bigl[a_j \cos(\omega j t) - b_j \sin(\omega j t)\bigr] \\ &= -\sum_{j=1}^J \omega^3 j^3 \bigl[a_j \cos(\omega j t) - b_j \sin(\omega j t)\bigr] \end{aligned} \]

직관: 사인·코사인 함수를 한 번 미분하면 주파수 \(\omega j\) 가 앞에 곱해지고, 세 번 미분하면 \((\omega j)^3\) 이 곱해진다.

Step 2: 연산자 \(L\) 적용

\[ \begin{aligned} L(x)(t) &= \omega^2 x^{(1)}(t) + x^{(3)}(t) \\ &= \sum_{j=1}^J \bigl[\omega^3 j - \omega^3 j^3\bigr] \bigl[a_j \cos(\omega j t) - b_j \sin(\omega j t)\bigr] \\ &= \sum_{j=1}^J \omega^3 j(1 - j^2) \bigl[a_j \cos(\omega j t) - b_j \sin(\omega j t)\bigr] \end{aligned} \]

직관: \(j = 1\) 일 때 \(1 - j^2 = 0\) 이므로, 기본 주파수 성분은 벌점을 전혀 받지 않는다. 이것이 조화 가속 연산자의 핵심 설계 의도이다 — 데이터의 기본 주기 패턴은 보존하고, 고주파 잡음만 벌점을 가한다.

Step 3: 적분 계산

\(\bigl[L(x)(t)\bigr]^2\) 를 적분할 때, 삼각함수의 직교성(orthogonality)을 활용한다:

\[ \int_0^T \cos(\omega j t) \cos(\omega k t) \, dt = \int_0^T \sin(\omega j t) \sin(\omega k t) \, dt = \frac{T}{2} \delta_{jk} \]

\[ \int_0^T \sin(\omega j t) \cos(\omega k t) \, dt = 0 \]

여기서 \(\delta_{jk}\) 는 크로네커 델타이다. 따라서 교차항은 모두 사라지고:

\[ \begin{aligned} \int_0^T \bigl[L(x)(t)\bigr]^2 \, dt &= \sum_{j=1}^J \omega^6 j^2 (1 - j^2)^2 \cdot \frac{T}{2} (a_j^2 + b_j^2) \\ &= \frac{T \omega^6}{2} \sum_{j=2}^J j^2 (j^2 - 1)^2 (a_j^2 + b_j^2) \end{aligned} \]

\(j = 1\) 항은 \((1 - 1)^2 = 0\) 이므로 합에서 제외된다. \((1-j^2)^2 = (j^2-1)^2\) 이므로 부호는 무관하다.

Step 4: \(T/2\) 와 \(\omega\) 관계 정리

\(\omega = 2\pi/T\) 이므로 \(T = 2\pi/\omega\) 이다. 따라서:

\[ \frac{T \omega^6}{2} = \frac{2\pi}{2\omega} \cdot \omega^6 = \pi \omega^5 \]

최종 결과:

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

1.3 해석

이 결과가 말하는 것:

\(j = 1\) (기본 주파수): 벌점 = 0 — 기본 연주기 패턴 보존
\(j = 2\): 벌점 계수 = \(4 \cdot 9 = 36\)
\(j = 3\): 벌점 계수 = \(9 \cdot 64 = 576\)
\(j = 10\): 벌점 계수 = \(100 \cdot 9801 = 980100\)

고주파로 갈수록 \(j^2(j^2-1)^2\) 이 폭발적으로 증가하므로, \(\lambda\) 가 작더라도 고주파 진폭 \(a_j, b_j\) 는 강하게 억제된다. 이것이 조화 가속 벌점이 주기적 데이터에 적합한 이유이다.

2 문제 2.2: FedYieldcurve 벌점 스무딩

2.1 문제

R 패키지 fds의 FedYieldcurve 데이터를 사용하여:

1982년 1월 이자율을 4개 B-spline 기저로 스무딩
6개 기저 + \(\lambda = 1\) + 2차 미분 벌점으로 재적합
다양한 \(\lambda\) 값 비교

2.2 풀이

library(fda)
library(fds)

# 데이터 추출 — 1982년 1월
data(FedYieldcurve)
maturities <- FedYieldcurve$x       # 만기(maturity)
jan82 <- FedYieldcurve$y[, "Jan-82"] # 1982년 1월 수익률

# (a) 4개 B-spline 기저로 단순 스무딩
basis4 <- create.bspline.basis(range(maturities), nbasis = 4)
fd4 <- smooth.basis(maturities, jan82, basis4)$fd

# (b) 6개 기저 + 벌점 스무딩 (lambda = 1, 2차 미분)
basis6 <- create.bspline.basis(range(maturities), nbasis = 6)
fdPar6 <- fdPar(basis6, Lfdobj = 2, lambda = 1)
fd6_pen <- smooth.basis(maturities, jan82, fdPar6)$fd

# 시각화
plot(maturities, jan82, pch = 19, xlab = "Maturity", ylab = "Yield (%)",
     main = "FedYieldcurve: Jan 1982")
lines(fd4, col = "blue", lwd = 2)
lines(fd6_pen, col = "red", lwd = 2)
legend("bottomright", c("Raw", "4 basis (no penalty)", "6 basis (lambda=1)"),
       pch = c(19, NA, NA), lty = c(NA, 1, 1), col = c("black", "blue", "red"))

(a) 해석: 4개 기저 함수는 기저 수가 적어 곡선의 자유도가 제한된다. 수익률 곡선의 전반적 형태(단기 고금리, 장기 저금리 — 역수익률 곡선)는 포착하지만, 중간 만기의 미세한 변곡점은 놓칠 수 있다.

(b) 해석: 6개 기저(데이터 점 수와 동일)를 사용하되 벌점을 부과하면, 이론적으로 데이터를 완벽히 보간할 수 있는 자유도를 가지면서도 \(\lambda = 1\) 이 과적합을 억제한다. 2차 미분 벌점은 “곡률이 작은” 매끄러운 곡선을 선호한다.

# (c) 여러 lambda 비교
lambdas <- c(0.001, 0.01, 0.1, 1, 10, 100)

plot(maturities, jan82, pch = 19, xlab = "Maturity", ylab = "Yield (%)")
colors <- rainbow(length(lambdas))

for (i in seq_along(lambdas)) {
  fdPar_i <- fdPar(basis6, Lfdobj = 2, lambda = lambdas[i])
  fd_i <- smooth.basis(maturities, jan82, fdPar_i)$fd
  lines(fd_i, col = colors[i], lwd = 2)
}

legend("bottomright", paste("lambda =", lambdas), col = colors, lwd = 2)

(c) 해석:

\(\lambda\)	동작	결과
매우 작음 (0.001)	벌점 무시 → 보간에 가까움	데이터 잡음까지 따라감
중간 (0.1~1)	적절한 균형	주요 패턴 보존, 잡음 제거
매우 큼 (100)	벌점 지배 → \(x''(t) \approx 0\)	직선에 수렴

1982년 1월 데이터는 점이 6개뿐이므로, \(\lambda \approx 0.1\) ~ \(1\) 범위가 가장 정보량이 풍부한 곡선을 생성한다. GCV를 사용하면 이 범위의 최적 \(\lambda\) 를 자동 선택할 수 있다.

3 문제 2.3: 상수 계수 미분 연산자

3.1 문제

계수 함수가 상수인 선형 미분 연산자:

\[ L(x)(t) = \sum_{m=0}^{M} \alpha_m x^{(m)}(t) \]

에 대해, 식 (2.1) 을 최소화하는 \(\hat{c}^\top = (\hat{c}_1, \ldots, \hat{c}_K)\) 의 표현식을 유도하라.

3.2 풀이

벌점 제곱합 (PSS):

\[ \mathrm{PSS}_\lambda(c_1, \ldots, c_K) = \sum_j \Bigl(y_j - \sum_{k=1}^K c_k B_k(t_j)\Bigr)^2 + \lambda \int_{T_1}^{T_2} \bigl[L(x_K)(t)\bigr]^2 \, dt \]

행렬 표기 도입:

\(\mathbf{c} = (c_1, \ldots, c_K)^\top\): 기저 계수 벡터
\(\boldsymbol{\Phi}\): \(J \times K\) 행렬, \(\Phi_{jk} = B_k(t_j)\) — 기저 함수의 관측점 평가 행렬
\(\mathbf{y} = (y_1, \ldots, y_J)^\top\): 관측값 벡터

그러면 적합 오차항은:

\[ \sum_j (y_j - x_K(t_j))^2 = \|\mathbf{y} - \boldsymbol{\Phi} \mathbf{c}\|^2 \]

벌점항의 행렬화:

\(L(x_K)(t) = \sum_{m=0}^M \alpha_m x_K^{(m)}(t)\) 이고, \(x_K^{(m)}(t) = \sum_k c_k B_k^{(m)}(t)\) 이므로:

\[ L(x_K)(t) = \sum_k c_k \Bigl[\sum_{m=0}^M \alpha_m B_k^{(m)}(t)\Bigr] \]

\(K \times K\) 벌점 행렬 \(\mathbf{R}\) 을 다음과 같이 정의한다:

\[ R_{k\ell} = \int_{T_1}^{T_2} \Bigl[\sum_{m=0}^M \alpha_m B_k^{(m)}(t)\Bigr] \Bigl[\sum_{m=0}^M \alpha_m B_\ell^{(m)}(t)\Bigr] \, dt \]

그러면 벌점항은:

\[ \int \bigl[L(x_K)(t)\bigr]^2 \, dt = \mathbf{c}^\top \mathbf{R} \, \mathbf{c} \]

최소화:

\[ \mathrm{PSS}_\lambda = \|\mathbf{y} - \boldsymbol{\Phi}\mathbf{c}\|^2 + \lambda \, \mathbf{c}^\top \mathbf{R} \, \mathbf{c} \]

\(\mathbf{c}\) 에 대해 미분하여 0으로 놓으면:

\[ -2 \boldsymbol{\Phi}^\top (\mathbf{y} - \boldsymbol{\Phi}\mathbf{c}) + 2\lambda \mathbf{R} \, \mathbf{c} = \mathbf{0} \]

정리하면:

\[ \boxed{\hat{\mathbf{c}} = (\boldsymbol{\Phi}^\top \boldsymbol{\Phi} + \lambda \mathbf{R})^{-1} \boldsymbol{\Phi}^\top \mathbf{y}} \]

3.3 해석

이 결과는 능선 회귀(ridge regression)와 구조적으로 동일하다:

능선 회귀: \(\hat{\boldsymbol{\beta}} = (\mathbf{X}^\top\mathbf{X} + \lambda \mathbf{I})^{-1} \mathbf{X}^\top \mathbf{y}\)
벌점 스무딩: \(\hat{\mathbf{c}} = (\boldsymbol{\Phi}^\top\boldsymbol{\Phi} + \lambda \mathbf{R})^{-1} \boldsymbol{\Phi}^\top \mathbf{y}\)

차이는 단위 행렬 \(\mathbf{I}\) 대신 벌점 행렬 \(\mathbf{R}\) 이 사용된다는 점이다. \(\mathbf{R}\) 은 미분 연산자의 구조를 반영하여 방향 선택적 축소(shrinkage)를 수행한다. 즉, 모든 계수를 균일하게 줄이는 것이 아니라, 고주파 성분의 계수를 더 강하게 줄인다.

4 문제 2.4: DTI 데이터의 등록 효과

4.1 문제

§1.5의 DTI 데이터에 대해:
(a) 100개 기저 + GCV로 벌점 스무딩
(b) 연속 등록 적용
(c) 등록 전후 FPCA 비교
(d) 등록의 필요성 판단

4.2 풀이

library(fda)
library(refund)

# DTI 데이터 로드 (refund 패키지)
data(DTI)
# FA 프로파일: tract positions × subjects
cca <- DTI$cca  # corpus callosum
tract <- seq(0, 1, length.out = ncol(cca))

# (a) 벌점 스무딩 (100 기저 + GCV)
basis100 <- create.bspline.basis(c(0, 1), nbasis = 100, norder = 6)
fdPar_gcv <- fdPar(basis100, Lfdobj = 4, lambda = 0)  # 초기값

# GCV 최적화
lambda_grid <- 10^seq(-10, -2, length.out = 20)
gcv_vals <- numeric(length(lambda_grid))

for (i in seq_along(lambda_grid)) {
  fdPar_i <- fdPar(basis100, Lfdobj = 4, lambda = lambda_grid[i])
  smooth_i <- smooth.basis(tract, t(cca), fdPar_i)
  gcv_vals[i] <- mean(smooth_i$gcv)
}

best_lambda <- lambda_grid[which.min(gcv_vals)]
fdPar_best <- fdPar(basis100, Lfdobj = 4, lambda = best_lambda)
dti_smooth <- smooth.basis(tract, t(cca), fdPar_best)

# 시각화
plot(dti_smooth$fd, xlab = "Tract Position", ylab = "FA",
     main = "DTI: Penalized Smoothing (GCV)")
mean_fd <- mean(dti_smooth$fd)
lines(mean_fd, lwd = 4, col = "black")

# (b) 연속 등록
regList <- register.fd(yfd = dti_smooth$fd)
dti_reg <- regList$regfd

plot(dti_reg, xlab = "Tract Position", ylab = "FA",
     main = "DTI: After Continuous Registration")
mean_reg <- mean(dti_reg)
lines(mean_reg, lwd = 4, col = "black")

# (c) FPCA 비교
pca_unreg <- pca.fd(dti_smooth$fd, nharm = 4)
pca_reg <- pca.fd(dti_reg, nharm = 4)

cat("== 비등록 FPCA ==\n")
cat("설명 분산 비율:", round(pca_unreg$varprop, 3), "\n")
cat("== 등록 후 FPCA ==\n")
cat("설명 분산 비율:", round(pca_reg$varprop, 3), "\n")

(a) 해석: 100개 B-spline 기저는 DTI FA 프로파일의 복잡한 형태를 충분히 포착할 수 있는 유연성을 제공한다. GCV가 최적 \(\lambda\) 를 선택하여, 측정 잡음은 제거하면서 백질 구조의 실제 FA 변동은 보존한다. 직접 스플라인 전개(Ch.1의 Figure 1.12)와 비교하면, 벌점 스무딩이 더 매끄러운 곡선을 생산한다.

(b) 해석: DTI 데이터에서 위상 변동은 “해부학적 랜드마크의 위치 차이”에 해당한다. 연속 등록은 각 피험자의 뇌량(corpus callosum) 내 구조적 특징점을 시간(tract position) 축을 따라 정렬하여, 개인 간 해부학적 차이를 보정한다. 등록 후 평균 함수가 더 선명한 봉우리를 보이는지 확인한다.

(c) 해석: 등록이 효과적이었다면:

첫 번째 FPC의 설명 분산 비율이 증가한다 (위상 변동이 제거되어 진폭 패턴에 집중)
더 적은 수의 FPC로 전체 변동을 설명할 수 있다
FPC의 형태가 더 해석 가능해진다 (위상 이동 패턴이 아닌 실제 구조 패턴 포착)

(d) 판단: DTI 데이터에서 등록의 필요성은 상황에 따라 다르다:

등록이 유용한 경우: 연구 목적이 “FA의 크기 차이”에 있을 때. 예: 질병군 vs 정상군의 FA 진폭 차이 탐지
등록이 불필요하거나 해로운 경우: 연구 목적이 “해부학적 구조의 위치 차이”에 있을 때. 등록이 바로 그 정보를 제거해 버리기 때문이다

DTI 데이터의 위상 변동은 대부분 개인의 뇌 크기·형태 차이에서 비롯된다. 이를 “잡음”으로 보고 제거할지, “신호”로 보고 분석 대상에 포함할지는 연구 질문에 달려 있다.

5 문제 2.5: 곡선 정렬의 위험성

5.1 문제

Bump 함수를 사용하여 곡선 정렬이 통계적 검정력을 해칠 수 있음을 보인다.

5.2 Bump 함수 정의

정의: Bump Function

중심 \(c_0\), 반지름 \(r_0\), 진폭 \(a_0\) 의 bump 함수:

\[ f(x) = \begin{cases} a_0 \exp\Bigl\{-\bigl(1 - \bigl(\frac{x - c_0}{r_0}\bigr)^2\bigr)^{-1}\Bigr\} & \text{if } |x - c_0| < r_0 \\ 0 & \text{otherwise} \end{cases} \]

이 함수는 무한히 미분 가능하면서도 지지(support)가 유한한 “부드러운 종 모양”이다. 직관적으로, 가우시안처럼 생겼지만 \(|x - c_0| \ge r_0\) 밖에서는 정확히 0이 된다.

5.3 (a) 시뮬레이션

library(fda)
library(MASS)  # mvrnorm

# Bump 함수 정의
bump <- function(x, c0, r0, a0) {
  result <- numeric(length(x))
  inside <- abs(x - c0) < r0
  z <- ((x[inside] - c0) / r0)^2
  result[inside] <- a0 * exp(-1 / (1 - z))
  result
}

# Matern 공분산 함수
matern_cov <- function(d, nu = 2.5, sigma = 1, rho = 0.1) {
  if (d == 0) return(sigma^2)
  x <- sqrt(2 * nu) * d / rho
  sigma^2 * (2^(1 - nu) / gamma(nu)) * x^nu * besselK(x, nu)
}

# 시뮬레이션 설정
n_grid <- 101
t_grid <- seq(0, 1, length.out = n_grid)
N <- 50  # 총 표본 크기

# 공분산 행렬 구성
Sigma <- matrix(0, n_grid, n_grid)
for (i in 1:n_grid) {
  for (j in 1:n_grid) {
    Sigma[i, j] <- matern_cov(abs(t_grid[i] - t_grid[j]))
  }
}

# 평균 함수: 그룹 1 (bump at 3/8), 그룹 2 (bump at 5/8)
mu1 <- bump(t_grid, c0 = 3/8, r0 = 1/4, a0 = 5)
mu2 <- bump(t_grid, c0 = 5/8, r0 = 1/4, a0 = 5)

# 함수형 표본 생성
set.seed(42)
curves <- matrix(0, N, n_grid)
group <- c(rep(0, N/2), rep(1, N/2))

for (i in 1:(N/2)) {
  curves[i, ] <- mu1 + mvrnorm(1, rep(0, n_grid), Sigma)
}
for (i in (N/2 + 1):N) {
  curves[i, ] <- mu2 + mvrnorm(1, rep(0, n_grid), Sigma)
}

# 시각화
matplot(t_grid, t(curves), type = "l", lty = 1,
        col = c(rep("blue", 25), rep("red", 25)),
        xlab = "t", ylab = "X(t)", main = "Simulated Bump Functions")
overall_mean <- colMeans(curves)
lines(t_grid, overall_mean, lwd = 4, col = "black")

핵심 관찰: 두 그룹의 곡선은 동일한 모양(bump)이지만, 위치가 다르다 — 그룹 1은 \(3/8\) 에, 그룹 2는 \(5/8\) 에 중심을 둔다. 전체 평균은 두 bump의 중간 (\(1/2\) 근처)에 넓고 낮은 봉우리를 보인다.

5.4 (b) 연속 등록 적용

# 함수형 객체 변환
basis <- create.bspline.basis(c(0, 1), nbasis = 30)
fdPar_obj <- fdPar(basis, Lfdobj = 2, lambda = 1e-4)
fd_obj <- smooth.basis(t_grid, t(curves), fdPar_obj)$fd

# 연속 등록
regList <- register.fd(yfd = fd_obj)
fd_reg <- regList$regfd

plot(fd_reg, xlab = "t", ylab = "X(t)",
     main = "After Continuous Registration")
mean_reg <- mean(fd_reg)
lines(mean_reg, lwd = 4, col = "black")

등록의 효과: 연속 등록은 모든 bump를 같은 위치로 이동시킨다. 그룹 1의 bump (\(3/8\))와 그룹 2의 bump (\(5/8\))가 모두 중간 지점으로 정렬된다. 결과적으로, 두 그룹 간의 유일한 차이(bump 위치)가 제거된다.

5.5 (c) FPCA + 회귀 분석

# 비등록 곡선: FPCA (1 PC)
pca_unreg <- pca.fd(fd_obj, nharm = 1)
scores_unreg <- pca_unreg$scores[, 1]

# 등록 곡선: FPCA (1 PC)
pca_reg <- pca.fd(fd_reg, nharm = 1)
scores_reg <- pca_reg$scores[, 1]

# 그룹 더미 변수에 대한 회귀
fit_unreg <- lm(scores_unreg ~ group)
fit_reg <- lm(scores_reg ~ group)

cat("== 비등록 ==\n")
cat("p-value:", summary(fit_unreg)$coefficients[2, 4], "\n")
cat("== 등록 후 ==\n")
cat("p-value:", summary(fit_reg)$coefficients[2, 4], "\n")

예상 결과 해석:

상태	p-value	해석
비등록	매우 작음 (유의)	첫 PC가 bump 위치 차이를 포착 → 그룹 차이 탐지
등록 후	큼 (비유의)	등록이 위치 차이를 제거 → 그룹 차이 소멸

직관: 비등록 상태에서 첫 번째 주성분은 “bump가 왼쪽에 있는가 오른쪽에 있는가”를 잡아낸다. 이 PC score에 그룹 더미를 회귀하면 당연히 유의하다. 그러나 등록 후에는 모든 bump가 같은 위치로 이동했으므로, PC가 위치 정보를 포착하지 못하고, 회귀의 유의성이 사라진다.

이것이 곡선 정렬의 위험이다: 위상 변동이 실제로 분석의 관심 대상일 때, 등록은 신호를 제거한다.

5.6 (d) 정렬이 해로운 실무 시나리오

시나리오: RT-PCR 증폭 곡선의 Ct 값 분석

RT-PCR에서 증폭 곡선의 모양은 모든 시료에서 유사한 시그모이드(S자)이다. 그러나 곡선이 역치를 넘는 시점(Ct 값)은 시료 내 초기 DNA 양에 비례한다:

높은 초기 DNA 양 → 빠른 증폭 → 작은 Ct
낮은 초기 DNA 양 → 느린 증폭 → 큰 Ct

이 상황에서 Ct 차이는 정확히 위상 변동이다. 만약 증폭 곡선을 정렬하면:

모든 곡선의 Ct가 같아진다
초기 DNA 양에 대한 정보가 완전히 사라진다
바이러스 정량, 유전자 발현 분석 등이 불가능해진다

교훈: 정렬은 도구이지 목적이 아니다. “위상 변동이 잡음인가 신호인가”를 반드시 먼저 판단해야 한다.

6 핵심 요약

문제	핵심 개념	실무 시사점
2.1	삼각함수 직교성 + 조화 가속 벌점	기본 주파수 보존, 고주파 억제
2.2	\(\lambda\) 선택의 편향-분산 트레이드오프	GCV로 자동 선택, 시각적 검증 필수
2.3	벌점 스무딩 = 방향 선택적 능선 회귀	\(\mathbf{R}\) 행렬이 축소 방향 결정
2.4	등록의 효과는 연구 질문에 의존	위상 vs 진폭 중 무엇이 관심인지 판단
2.5	정렬이 신호를 제거할 수 있음	위상 변동이 신호이면 정렬 금지

7 관련 주제

선행 지식

후속 주제

관련 개념