Kwangmin Kim - FDA 1.1~1.2 — 기저 전개와 표본 평균·공분산

1 기저 전개 (Basis Expansion)

1.1 왜 기저 전개가 필요한가

함수형 데이터의 원시 형태는 이산 관측값의 집합이다:

\[ x_n(t_{j,n}) \in \mathbb{R}, \quad t_{j,n} \in [T_1, T_2], \quad n = 1, \ldots, N, \quad j = 1, \ldots, J_n \]

$N$ 개의 곡선이 공통 구간 $[T_1, T_2]$ 위에서 관측되지만, 각 곡선 $x_n$ 의 관측 시점 $t_{j,n}$ 은 곡선마다 다를 수 있고, 관측점 사이의 값은 알 수 없다. 이 상태 그대로는 곡선 간 비교도, 미분도, 적분도 불가능하다.

기저 전개(basis expansion)는 이 이산 데이터를 연속 함수 객체로 변환하는 핵심 기법이다:

\[ X_n(t) \approx \sum_{m=1}^{M} c_{nm} B_m(t), \quad 1 \leq n \leq N \]

여기서 $B_m$ 은 미리 정해진 기저 함수(basis function), $c_{nm}$ 은 곡선 $n$ 의 $m$ 번째 전개 계수(expansion coefficient)이다. $M$ 은 기저 함수의 개수로, 보통 관측점 수 $J_n$ 보다 작다.

이 전개가 하는 일을 비유하면, “10,000개 픽셀로 된 사진을 25개 핵심 패턴의 가중합으로 압축하는 것”과 유사하다. 원래 사진의 세부 잡음은 버리면서도 전체 형태는 보존한다.

1.2 기저 전개의 세 가지 역할

기저 전개는 단순한 커브 피팅(curve fitting)이 아니라, FDA의 전 과정을 가능하게 하는 기반 인프라이다:

역할 1: 차원 축소

자기장 관측소의 하루 데이터가 17,280개 점(5초 간격)이라면, $M = 25$ 개의 B-spline 기저로 전개하면 25개 계수만으로 하루의 곡선을 표현할 수 있다. 이는 692:1 압축이다. 고빈도 금융 데이터(분당 1점, 하루 390점)도 마찬가지이다.

각 곡선 $x_n$ 은 $M$ 차원 벡터 $\mathbf{c}_n = [c_{n1}, c_{n2}, \ldots, c_{nM}]^T$ 로 요약된다. 이 계수 벡터에 기존의 다변량 통계 기법(PCA, 회귀, 군집 분석 등)을 적용할 수 있다는 실용적 이점이 생긴다. 실제로 많은 FDA 조작은 기저 계수에 대한 행렬 연산으로 귀결된다 — 이것이 FDA를 계산적으로 tractable하게 만드는 핵심 메커니즘이다.

역할 2: 공통 도메인 확보

관측 시점 $t_{j,n}$ 이 곡선마다 다른 경우가 흔하다. 예를 들어, 환자 A는 1월·4월·9월에 혈액 검사를 받고, 환자 B는 2월·6월·11월에 검사를 받았다면, 두 곡선의 관측 시점이 다르다. 이 상태에서는 두 곡선의 3월 값을 직접 비교할 수 없다. 기저 전개를 거치면 두 곡선 모두 연속 함수가 되어, 임의의 시점 $t$ 에서 값을 계산하고 비교할 수 있다. 이것이 다변량 분석의 결측치 문제를 FDA가 우회하는 방법이다.

역할 3: 예비 평활화(Preliminary Smoothing)

기저 함수가 매끄러우면(B-spline, Fourier 등) 전개 결과도 매끄럽다. 원시 데이터에 포함된 측정 잡음이 자연스럽게 제거되는 효과가 있다. 물론 더 정교한 평활화(벌점 평활 등)는 이후에 별도로 적용할 수 있다 (Ch.2에서 다룬다).

1.3 B-spline 기저

B-spline(Basis spline)은 FDA에서 가장 범용적으로 사용되는 기저 시스템이다. 각 B-spline 기저 함수는 구간의 일부에서만 0이 아닌 값을 가진다 — 이를 국소 지지(local support)라 한다.

B-spline의 핵심 성질

국소 지지: 각 기저 함수가 전체 구간이 아닌 인접 매듭(knot) 사이에서만 활성화된다.
비음성성: 모든 $t$ 에서 $B_m(t) \geq 0$ 이다.
단위 분할: 모든 $t$ 에서 $\sum_m B_m(t) = 1$ 이다.
차수(order) $k$: $k-2$ 차 연속 도함수를 보장한다. 기본값 $k = 4$ (3차 스플라인)이면 2차 도함수까지 연속이다.

국소 지지가 왜 중요한가? 곡선의 한 구간에서 날카로운 변화가 발생해도 다른 구간의 적합에 영향을 주지 않는다. 이는 Fourier 기저와의 근본적인 차이점이다 — Fourier 기저 함수는 전체 구간에 걸쳐 진동하므로, 한 지점의 이상이 곡선 전체에 파급된다.

R에서 B-spline 기저를 생성하는 코드는 다음과 같다:

library(fda)
spline.basis <- create.bspline.basis(rangeval = c(0, 10), nbasis = 5)
plot(spline.basis, lty = 1, lwd = 2)

create.bspline.basis() 의 주요 인수:

인수	의미	기본값
`rangeval`	정의 구간 $[T_1, T_2]$	필수
`nbasis`	기저 함수 수 $M$	필수
`norder`	차수 $k$ (스플라인 차수 + 1)	4 (3차 스플라인)
`breaks`	내부 매듭 위치	등간격 자동 배치

nbasis 와 norder, breaks 사이에는 관계 $M = \text{norder} + \text{nbreaks} - 2$ 가 성립한다. 매듭을 많이 넣으면 유연성이 증가하지만, 과적합(overfitting) 위험이 커진다. 이 trade-off는 벌점 평활(Ch.2)에서 체계적으로 다룬다.

1.4 Fourier 기저

Fourier 기저는 상수 함수와 주기가 감소하는(주파수가 증가하는) 사인·코사인 쌍으로 구성된다. $[0, T]$ 구간에서 $2K+1$ 개의 Fourier 기저 함수는:

\[ 1, \quad \sin\left(\frac{2\pi t}{T}\right), \quad \cos\left(\frac{2\pi t}{T}\right), \quad \sin\left(\frac{4\pi t}{T}\right), \quad \cos\left(\frac{4\pi t}{T}\right), \quad \ldots \]

첫 번째 함수는 상수(평균 수준), 이후 사인/코사인 쌍이 점차 높은 주파수의 변동을 포착한다. Fourier 기저의 특징은 항상 홀수 개($2K+1$)를 사용한다는 것이다 — 상수 함수 1개와 사인/코사인 쌍 $K$ 개이다.

fourier.basis <- create.fourier.basis(rangeval = c(0, 10), nbasis = 5)
plot(fourier.basis, lty = 1, lwd = 2)

Fourier 기저는 주기적 데이터에 이상적이다: 연간 기온 변화, 일별 전력 소비 패턴, 심전도(ECG) 신호 등. 그러나 한 가지 중요한 제약이 있다:

Fourier 기저의 제약

Fourier 기저는 구간의 양 끝에서 함수값이 비슷한 데이터에만 적합하다. 사인·코사인은 본질적으로 주기 함수이므로, $x(T_1) \neq x(T_2)$ 인 데이터에 Fourier 기저를 적용하면 양 끝에서 인위적인 진동이 발생할 수 있다 (Gibbs 현상).

이 때문에 Fourier 기저는 비주기적 데이터(주가, 성장 곡선 등)에는 적합하지 않다. BOA 주식 수익률 데이터의 경우, 매일 0에서 시작하지만 종가가 0이 아니므로 양 끝 값이 다르다 — 이런 데이터에는 B-spline이 더 적합하다.

1.5 B-spline vs Fourier: 선택 기준

특성	B-spline	Fourier
지지(support)	국소적 — 인접 매듭 범위	전역적 — 전체 구간
주기성	비주기적 데이터에 적합	주기적 데이터에 최적
경계 행동	양 끝 독립적 처리	$x(T_1) \approx x(T_2)$ 가정
미분 가능 횟수	차수 $k$ 에 의해 제한 ($k-2$ 회)	무한 미분 가능
국소 변화 반응	국소적으로만 영향	전체 곡선에 파급
매듭 배치	유연 — 변화가 많은 구간에 밀집 가능	해당 없음 (전역 주파수)

실무적 경험 규칙: 데이터의 주기성 여부가 확실하지 않으면 B-spline을 기본으로 선택한다. B-spline은 주기적 데이터에도 합리적으로 작동하지만(매듭을 충분히 넣으면), Fourier는 비주기적 데이터에서 경계 문제를 일으킬 수 있기 때문이다.

1.6 예시: 위너 과정의 B-spline 전개

위너 과정(Wiener process, 브라운 운동)은 랜덤 워크의 연속 시간 극한이다. 구간 $[0, K]$ 에서 랜덤 워크를 다음과 같이 정의한다:

\[ S_i = \frac{1}{\sqrt{K}} \sum_{k=1}^{i} N_k, \quad N_k \stackrel{\text{iid}}{\sim} N(0, 1), \quad 1 \leq k \leq K \]

각 단계에서 표준정규 변량 $N_k$ 를 누적하고, $\sqrt{K}$ 로 스케일링하여 분산을 정규화한다. $K = 10{,}000$ 일 때, 이 랜덤 워크를 함수 $X(t_i) = S_i$ ($t_i = i$)로 간주하면 10,000개의 점으로 이루어진 하나의 “곡선”이 된다.

이 곡선을 25개의 B-spline 기저 함수로 전개하면:

library(fda)
Wiener <- cumsum(rnorm(10000)) / 100
plot.ts(Wiener, xlab = "", ylab = "")
B25.basis <- create.bspline.basis(rangeval = c(0, 10000), nbasis = 25)
Wiener.fd <- smooth.basis(y = Wiener, fdParobj = B25.basis)
lines(Wiener.fd, lwd = 3)

smooth.basis() 는 원시 데이터 y 와 기저 객체 fdParobj 를 받아 함수 데이터 객체(fd object)를 반환한다. 반환값은 전개 계수 $c_{nm}$ 과 기저 정보를 모두 포함하는 리스트이다. 이후의 모든 FDA 연산(평균, 분산, PCA, 회귀 등)은 이 fd 객체를 입력으로 받는다.

그래프에서 관찰할 수 있는 핵심:

회색의 원래 랜덤 워크는 10,000개의 점으로 이루어져 있고 잡음이 크다
굵은 선으로 그린 B-spline 전개는 25개 계수만으로 전체 추세를 충실히 포착한다
세밀한 등락(고주파 잡음)은 제거되었지만, 전체적인 상승/하강 패턴은 보존된다
이것이 400:1 차원 축소($10{,}000 \to 25$)이면서 동시에 예비 평활화인 것이다

1.7 기저 수 $M$ 의 선택

$M$ 이 너무 작으면 곡선의 중요한 특징을 놓치고(underfitting), 너무 크면 잡음까지 적합한다(overfitting). 이 trade-off를 다루는 체계적 방법이 벌점 평활(penalized smoothing)이며, Ch.2에서 상세히 다룬다. 이 장에서는 일반화 교차 검증(GCV)을 통해 $M$ 또는 스무딩 모수 $\lambda$ 를 선택한다는 정도만 언급한다.

경험적으로, B-spline의 경우 관측점 수의 제곱근 정도($M \approx \sqrt{J}$)가 합리적 출발점이다. 위의 예시에서 $J = 10{,}000$ 이고 $M = 25 \approx \sqrt{625}$ 를 사용했는데, 이는 데이터 복잡도에 비해 다소 보수적인 선택이었다. 더 많은 기저를 쓰면 세밀한 변동을 포착할 수 있지만, 랜덤 워크의 경우 그 세밀한 변동 자체가 잡음이므로 25개가 적절하다.

1.8 추가적 평활 기법에 대한 언급

기저 전개 외에도, Ramsay et al. (2009)은 더 정교한 함수 객체 구성 방법을 다룬다:

벌점 평활(roughness penalty smoothing): 기저 전개에 $\int [x^{(k)}(t)]^2 dt$ 형태의 거칠기 벌점을 추가하여 과적합을 방지한다 (Ch.2의 주제)
단조 평활(monotone smoothing): 성장 곡선처럼 값이 감소할 수 없는 데이터에 사용한다. 기저 전개는 원칙적으로 $X_n(t)$ 의 단조성을 보장하지 않는다 — 실제로 성장 데이터에 단순 B-spline 전개를 적용하면 작은 구간에서 키가 줄어드는 비물리적 결과가 나올 수 있다
양수 평활(positivity-preserving smoothing): 농도, 확률 등 음수가 될 수 없는 양에 사용한다

이러한 제약 평활은 R의 fda 패키지에서 customized 함수로 지원된다.

2 표본 평균과 공분산

기저 전개를 통해 원시 데이터를 함수 객체로 변환한 뒤, 함수형 데이터의 가장 기본적인 요약 통계량을 계산할 수 있다.

2.1 점별 평균 함수 (Pointwise Sample Mean)

$N$ 개의 함수 관측 $X_1(t), X_2(t), \ldots, X_N(t)$ 에 대해, 표본 평균 함수는 각 시점 $t$ 에서의 산술 평균이다:

\[ \bar{X}_N(t) = \frac{1}{N} \sum_{n=1}^{N} X_n(t) \]

이것은 스칼라 표본 평균 $\bar{x} = \frac{1}{N}\sum x_i$ 의 함수 공간 확장이다. 스칼라 평균이 하나의 숫자인 것처럼, 함수 평균은 하나의 곡선이다.

직관적으로, $\bar{X}_N(t)$ 는 “전형적인 곡선의 형태”를 나타낸다. 100명의 환자의 식후 혈당 곡선을 평균하면, 식후 30분에 급상승 $\to$ 2시간 후 서서히 감소라는 전형적 패턴이 드러난다. 개별 환자의 곡선은 이 전형적 패턴 주변에서 변동한다.

2.2 점별 표준편차 함수 (Pointwise Sample SD)

\[ \text{SD}_X(t) = \left\{ \frac{1}{N-1} \sum_{n=1}^{N} \left( X_n(t) - \bar{X}_N(t) \right)^2 \right\}^{1/2} \]

$\text{SD}_X(t)$ 는 시점 $t$ 에서 곡선들이 평균 곡선으로부터 얼마나 흩어져 있는지를 보여준다. 이 함수가 시점마다 다른 값을 가진다는 것이 핵심이다.

예를 들어, 브라운 운동 50개의 표본에서:

$t$ 가 작을 때(시작점 근처): 모든 곡선이 0 근처이므로 $\text{SD}_X(t)$ 가 작다
$t$ 가 클 때(구간 끝 근처): 곡선들이 각자의 방향으로 크게 벌어져 있으므로 $\text{SD}_X(t)$ 가 크다
즉 $\text{SD}_X(t)$ 는 $t$ 에 대해 증가하는 함수이다 — 이는 브라운 운동의 분산이 $\text{Var}(W(t)) = t$ 로 시간에 비례하여 증가한다는 이론적 결과와 일치한다

이처럼 표준편차 함수는 “곡선 간 변동이 어느 구간에서 크고 어디서 작은지”를 시각적으로 파악하게 해준다. 스칼라 통계에서는 하나의 숫자(SD)만 있지만, FDA에서는 변동의 구간별 패턴을 볼 수 있다 — 이것이 FDA가 정보를 더 많이 추출하는 한 가지 방식이다.

2.3 예시: 50개 랜덤 워크의 평균과 SD

N <- 50
W.mat <- matrix(0, ncol = N, nrow = 10000)
for (n in 1:N) { W.mat[, n] <- cumsum(rnorm(10000)) / 100 }
B25.basis <- create.bspline.basis(rangeval = c(0, 10000), nbasis = 25)
W.fd <- smooth.basis(y = W.mat, fdParobj = B25.basis)
plot(W.fd, ylab = "", xlab = "", col = "gray", lty = 1)
W.mean <- mean(W.fd$fd)
W.sd <- std.fd(W.fd$fd)
lines(W.sd, lwd = 3)          # SD: 굵은 실선
lines(W.mean, lty = 2, lwd = 3)  # 평균: 굵은 점선

코드의 흐름:

matrix(0, ncol=N, nrow=10000): $50 \times 10{,}000$ 데이터 행렬 준비
cumsum(rnorm(10000)) / 100: 각 열에 랜덤 워크 생성
smooth.basis(): 50개 랜덤 워크를 일괄적으로 함수 객체로 변환
W.fd$fd: smooth.basis 의 반환값은 리스트이며, $fd 로 함수 값을 추출한다
mean(), std.fd(): 함수 객체에 대한 점별 평균과 SD를 계산한다

그래프에서 회색 곡선 50개가 중앙의 점선(평균 $\approx 0$) 주변으로 퍼져 있고, 굵은 실선(SD)이 $t$ 가 증가할수록 커지는 패턴을 확인할 수 있다.

2.4 표본 공분산 함수 (Sample Covariance Function)

점별 SD는 각 시점에서의 변동 크기만 알려줄 뿐, 서로 다른 시점 간의 관계는 담지 못한다. “시점 $t$ 에서 평균보다 높은 곡선이 시점 $s$ 에서도 높은가?” — 이 질문에 답하려면 공분산 함수가 필요하다:

\[ \hat{c}(t, s) = \frac{1}{N-1} \sum_{n=1}^{N} \left( X_n(t) - \bar{X}_N(t) \right) \left( X_n(s) - \bar{X}_N(s) \right) \]

공분산 행렬에서 공분산 함수로

다변량 분석에서 $p$ 차원 벡터의 공분산 행렬 $\hat{\Sigma}$ 는 $p \times p$ 이다. FDA에서는 벡터가 함수로 확장되므로, 공분산도 행렬이 아닌 $[T_1, T_2] \times [T_1, T_2]$ 위의 이변량 함수로 확장된다.

다변량	FDA
$\hat{\Sigma}_{jk} = \frac{1}{N-1}\sum_{n}(x_{nj} - \bar{x}_j)(x_{nk} - \bar{x}_k)$	$\hat{c}(t,s) = \frac{1}{N-1}\sum_{n}(X_n(t) - \bar{X}_N(t))(X_n(s) - \bar{X}_N(s))$
인덱스: 이산 $(j, k) \in \{1, \ldots, p\}^2$	인덱스: 연속 $(t, s) \in [T_1, T_2]^2$
시각화: 행렬 히트맵	시각화: 3D 표면(perspective plot) 또는 등고선도(contour plot)

2.5 공분산 함수의 해석

$\hat{c}(t, s)$ 의 해석은 스칼라 공분산과 동일하다:

$\hat{c}(t, s) > 0$: 시점 $t$ 에서 평균보다 높은 곡선은 시점 $s$ 에서도 평균보다 높은 경향이 있다 (양의 공변동)
$\hat{c}(t, s) < 0$: 시점 $t$ 에서 높은 곡선이 시점 $s$ 에서는 낮은 경향이 있다 (음의 공변동)
$\hat{c}(t, s) = 0$: 두 시점에서의 값이 (선형적으로) 무관하다
$\hat{c}(t, t) = \text{Var}(X(t))$: 대각선은 각 시점의 분산이다. 따라서 점별 SD 함수의 제곱은 공분산 함수의 대각선 값이다: $\text{SD}_X(t)^2 = \hat{c}(t, t)$

공분산 함수는 FDA에서 단순한 기술 통계를 넘어 핵심적 역할을 한다: 함수 주성분 분석(FPCA)은 공분산 함수의 고유분해를 통해 수행되고, 함수 회귀의 추론도 공분산 구조에 의존한다. 따라서 공분산 함수의 정확한 추정은 이후 모든 분석의 기초이다.

2.6 시각화: 3D 표면과 등고선도

공분산 함수 $\hat{c}(t, s)$ 는 이변량 함수이므로, 시각화에는 두 가지 방법이 있다:

Perspective plot (3D 표면): $(t, s)$ 평면 위에 $\hat{c}(t,s)$ 를 높이로 표현한다. 전체 구조를 직관적으로 파악하는 데 유용하다.
Contour plot (등고선도): 같은 공분산 값을 가진 $(t, s)$ 점들을 등고선으로 연결한다. 정밀한 비교와 비등방성 파악에 유용하다.

W.cov <- var.fd(W.fd$fd)
grid <- (1:100) * 100
W.cov.mat <- eval.bifd(grid, grid, W.cov)
persp(grid, grid, W.cov.mat, xlab = "s", ylab = "t", zlab = "c(s,t)")
contour(grid, grid, W.cov.mat, lwd = 2)

코드의 핵심:

var.fd(): 함수 데이터 객체의 공분산 함수를 계산한다. 결과는 bifd (이변량 함수 데이터) 객체이다
eval.bifd(): 이변량 함수 객체를 지정한 격자점 (grid, grid) 에서 평가하여 행렬로 반환한다
이 행렬이 시각화의 입력이 된다

2.7 브라운 운동의 이론적 공분산

50개 랜덤 워크(브라운 운동의 이산 근사)로부터 추정한 공분산 함수는 이론적 결과와 비교할 수 있다. 브라운 운동 $W(t)$ 의 모공분산 함수는 닫힌 형태로 알려져 있다:

\[ c(t, s) = \min(t, s) \]

이 결과의 직관적 이해: 브라운 운동은 독립 증분(independent increment)을 가지므로, 시점 $t < s$ 에서의 공분산은 두 시점이 “공유하는 이력”의 길이에 의해 결정된다. $W(t)$ 는 $[0, t]$ 구간의 증분의 합이고, $W(s)$ 는 $[0, s]$ 구간의 증분의 합이다. 두 값이 공유하는 구간은 $[0, t]$ (더 이른 시점까지)이므로, 공분산은 $\text{Cov}(W(t), W(s)) = t = \min(t, s)$ 이다.

등고선도에서 이 이론적 구조가 명확히 드러난다:

등고선이 대각선을 중심으로 삼각형 모양을 이룬다
대각선($t = s$)에서의 값 $c(t, t) = t$ 는 분산, 즉 $t$ 에 비례하여 증가한다
대각선에서 멀어질수록(시점 간 거리가 멀수록) 공분산이 감소한다 — 다만 $\min(t,s)$ 이므로 0으로 줄지는 않는다

추정된 공분산 표면이 이 이론적 $\min(t, s)$ 패턴과 유사하다는 것은, 기저 전개 $\to$ 공분산 추정이라는 FDA 파이프라인이 데이터의 진짜 구조를 잘 포착하고 있음을 보여준다.

2.8 점별 신뢰 구간

평균 함수에 대한 점별 95% 신뢰 구간은 다음과 같이 구성된다:

\[ \bar{X}_N(t) \pm 2 \cdot \frac{\text{SD}_X(t)}{\sqrt{N}} \]

이 구간은 각 시점 $t$ 에서 독립적으로 구성된다 — 따라서 “점별(pointwise)”이라 부른다. 이것은 “모든 $t$ 에서 동시에” 평균을 커버하는 동시 신뢰 대역(simultaneous confidence band)보다 좁다 (동시 신뢰 대역은 Ch.12에서 다룬다).

R에서의 구현은 기저 계수에 대한 산술로 이루어진다:

muhat <- mean.fd(log_BOA_f)
sdhat <- sd.fd(log_BOA_f)
SE_hat_U <- fd(basisobj = bspline_basis)
SE_hat_L <- fd(basisobj = bspline_basis)
SE_hat_U$coefs <- 2 * sdhat$coefs / sqrt(N) + muhat$coefs
SE_hat_L$coefs <- -2 * sdhat$coefs / sqrt(N) + muhat$coefs

여기서 핵심은 $coefs — 기저 전개의 계수 벡터 — 에 대한 스칼라 연산으로 함수 수준의 조작(평균 $\pm$ 2SE)을 수행한다는 것이다. 기저 전개가 모든 FDA 연산을 계수 벡터 연산으로 환원한다는 원리의 구체적 예이다.

3 코드 예시: Python 구현

3.1 Step 1: 순수 Python으로 기저 전개 원리 확인

import numpy as np
from scipy.interpolate import BSpline

np.random.seed(42)

K = 10000
M = 25
t = np.linspace(0, 1, K)

walk = np.cumsum(np.random.randn(K)) / np.sqrt(K)

knots_interior = np.linspace(0, 1, M - 2)
order = 4  # cubic B-spline
knots_augmented = np.concatenate([
    np.repeat(0.0, order),
    knots_interior[1:-1],
    np.repeat(1.0, order)
])
n_basis = len(knots_augmented) - order

B = np.zeros((K, n_basis))
for i in range(n_basis):
    coeffs = np.zeros(n_basis)
    coeffs[i] = 1.0
    spl = BSpline(knots_augmented, coeffs, order - 1, extrapolate=False)
    B[:, i] = spl(t)

c_hat = np.linalg.lstsq(B, walk, rcond=None)[0]
X_smooth = B @ c_hat

print(f"원시 데이터 점 수: {K}")
print(f"기저 수 (계수 수): {n_basis}")
print(f"압축 비율: {K / n_basis:.0f}:1")
print(f"잔차 제곱합: {np.sum((walk - X_smooth)**2):.4f}")

최소제곱법 np.linalg.lstsq(B, walk) 은 $\|\mathbf{y} - \mathbf{B}\mathbf{c}\|^2$ 를 최소화하는 계수 벡터 $\mathbf{c}$ 를 구한다. 이것이 기저 전개의 실제 연산이다 — 기저 행렬 $\mathbf{B}$ ($K \times M$)와 데이터 벡터 $\mathbf{y}$ ($K \times 1$)에 대한 선형 회귀이다.

3.2 Step 2: 평균·공분산 계산

import numpy as np
from scipy.interpolate import BSpline

np.random.seed(42)

N, K, M = 50, 10000, 25
t = np.linspace(0, 1, K)

walks = np.cumsum(np.random.randn(N, K), axis=1) / np.sqrt(K)

knots_interior = np.linspace(0, 1, M - 2)
order = 4
knots_augmented = np.concatenate([
    np.repeat(0.0, order),
    knots_interior[1:-1],
    np.repeat(1.0, order)
])
n_basis = len(knots_augmented) - order

B = np.zeros((K, n_basis))
for i in range(n_basis):
    coeffs = np.zeros(n_basis)
    coeffs[i] = 1.0
    spl = BSpline(knots_augmented, coeffs, order - 1, extrapolate=False)
    B[:, i] = spl(t)

C_all = np.linalg.lstsq(B, walks.T, rcond=None)[0].T  # (N, n_basis)
smoothed = C_all @ B.T  # (N, K): 평활된 곡선 행렬

mean_func = smoothed.mean(axis=0)
sd_func = smoothed.std(axis=0, ddof=1)

centered = smoothed - mean_func[np.newaxis, :]
cov_matrix = centered.T @ centered / (N - 1)  # (K, K)

print(f"평균 함수 범위: [{mean_func.min():.4f}, {mean_func.max():.4f}]")
print(f"SD 함수 범위: [{sd_func.min():.4f}, {sd_func.max():.4f}]")
print(f"공분산 행렬 shape: {cov_matrix.shape}")

t_mid = K // 2
print(f"\nc(t_mid, t_mid) = Var(X(t_mid)) = {cov_matrix[t_mid, t_mid]:.4f}")
print(f"SD(t_mid)^2 = {sd_func[t_mid]**2:.4f}")
print(f"이론값 Var(W(0.5)) = 0.5 = {0.5:.4f}")

마지막 세 줄이 핵심 검증이다: $\hat{c}(t, t)$ 는 해당 시점의 분산이며, 브라운 운동의 이론적 분산 $\text{Var}(W(t)) = t$ 와 비교할 수 있다. $t = 0.5$ 에서 추정 분산이 0.5 근처인지 확인한다.

3.3 Step 3: scikit-fda 구현 (실무 활용)

import skfda
from skfda.representation.basis import BSplineBasis
import numpy as np

np.random.seed(42)
N, K = 50, 10000
walks = np.cumsum(np.random.randn(N, K), axis=1) / np.sqrt(K)
t_grid = np.linspace(0, 1, K)

fd = skfda.FDataGrid(data_matrix=walks, grid_points=t_grid)
basis = BSplineBasis(n_basis=25, domain_range=(0, 1))
fd_basis = fd.to_basis(basis)

mean_fd = fd_basis.mean()
print("평균 함수 (기저 계수):", mean_fd.coefficients.flatten()[:5], "...")

var_fd = fd.var()
cov_matrix = np.cov(fd.data_matrix.squeeze(), rowvar=True)
print(f"공분산 행렬 shape: {cov_matrix.shape}")

skfda.FDataGrid 는 격자 관측 함수 데이터를, to_basis() 는 기저 표현으로의 변환을 수행한다. R의 smooth.basis() 와 동일한 역할이다.

4 관련 주제

선행 지식

FDA 시리즈

후속 주제

FDA 1.3 — 주성분 함수 (EFPC)
FDA 1.4~1.5 — BOA 주식 수익률·DTI 응용
FDA Ch.2 — 미분, 벌점 스무딩, 곡선 정렬