FDA 7.6 — Chapter 7 연습문제 풀이

Bias-Variance 분해 검증·MSE 최적화·Local polynomial 닫힌 해·RKHS Representer Theorem·BLUP 및 Taylor 전개의 도함수 추정

Kokoszka & Reimherr (2017) Ch.7 의 연습문제 17 개를 상세 풀이한다. Problem 7.1~7.4 는 sparse FDA 의 점근 분석 (Bias 식 7.2, Taylor 전개 7.3, 분산 식 7.5, MSE 최적 h 와 M = N^δ 분석), Problem 7.5~7.6 은 Local polynomial 의 닫힌 해 (7.9) 와 RKHS 노름 (7.10) 검증, Problem 7.7~7.10 은 R 코드 (RKHS, gam, CATT 공분산), Problem 7.11~7.13 은 RKHS 의 이론 (Representer Theorem, 최소화 형태, Canadian Weather 비교), Problem 7.14~7.17 은 양정치 보정·BLUP·Taylor 전개로 도함수 추정·Gaussian/Exponential 핵의 매끄러움 — 모두 직관·수식 유도·R 코드 포함.

Statistics
Functional Data Analysis
저자

Kwangmin Kim

공개

2026년 05월 07일

1 이 포스트의 위치

./7-0-sparse-fda-overview.qmd 부터 ./7-3-sparse-regression.qmd 까지 Ch.7 의 본문 5 개 절을 다뤘다. 이 포스트는 마지막 §7.6 의 17 개 연습문제를 풀이한다.

Ch.7 본문 (§7.1~§7.5)
    ↓
Ch.7 연습문제 (§7.6, 17 문제)
    ├ 7.1~7.4: 점근 분석 (Bias·Variance·MSE 최적·M=N^δ 시나리오)
    ├ 7.5~7.6: Local polynomial 닫힌 해 + RKHS 노름
    ├ 7.7~7.10: R 코드 (RKHS·gam·CATT·미국 기상 데이터)
    ├ 7.11~7.13: RKHS 이론 (Representer Theorem·최소화 형태·Canadian Weather)
    └ 7.14~7.17: 보정·BLUP·도함수 추정·핵 매끄러움

핵심 메시지: 연습문제는 7.1 의 점근 분석 결과 (M ~ N^{1/4} 임계값) 와 7.2 의 RKHS 도구 (Representer Theorem, 핵 선택의 매끄러움 효과) 의 정확한 검증. 특히 Problem 7.4 의 MSE 최적화와 Problem 7.16 의 Taylor 전개 도함수 추정이 sparse FDA 의 핵심 통찰이다.


2 Problem 7.1: 식 (7.2) 의 정당화

2.1 문제

\(N = \sum_n M_n\). 식 (7.2) \(E[\widehat{\mu}_h(t) \mid \{t_{nm}\}] \approx (2h)^{-1} \int_{t-h}^{t+h} \mu(x) \, dx\) 를 정당화하시오. \(h\) 에 어떤 조건이 필요한가?

2.2 직관

\(\widehat{\mu}_h(t)\)\(t\) 주위 데이터의 가중 평균이므로, 그 기댓값은 \(t\) 주위 평균 함수의 가중 평균. 박스 커널이면 단순한 균일 평균.

2.3 풀이

박스 커널 가정 \(K(x) = \mathbb{1}_{|x| \leq 1}\). Nadarya-Watson 추정량:

\[ \widehat{\mu}_h(t) = \frac{\sum_n \sum_m K\left(\frac{t - t_{nm}}{h}\right) Y_{nm}}{\sum_n \sum_m K\left(\frac{t - t_{nm}}{h}\right)}. \]

조건부 기댓값 (시점 \(\{t_{nm}\}\) 고정):

\[ E[\widehat{\mu}_h(t) \mid \{t_{nm}\}] = \frac{\sum_n \sum_m K\left(\frac{t - t_{nm}}{h}\right) E[Y_{nm} \mid t_{nm}]}{\sum_n \sum_m K\left(\frac{t - t_{nm}}{h}\right)} = \frac{\sum K \cdot \mu(t_{nm})}{\sum K}. \]

분자/분모 모두 표본 평균. \(t_{nm}\)\(U(0, 1)\) iid 이므로 (Example 7.1.2 가정), 큰 \(N = \sum_n M_n\) 의 한계에서:

  • 분모 \(\approx 2 N h\) (커널 안의 점 수, 박스 폭 \(2h\)).
  • 분자 \(\approx N \int_{t-h}^{t+h} \mu(x) \, dx\) (균일 측도 적분).

따라서:

\[ E[\widehat{\mu}_h(t) \mid \{t_{nm}\}] \approx \frac{N \int_{t-h}^{t+h} \mu(x) \, dx}{2 N h} = (2h)^{-1} \int_{t-h}^{t+h} \mu(x) \, dx. \]

2.4 조건

  • \(h \to 0\) (커널이 \(t\) 주위로 집중되어야 의미 있음).
  • \(Nh \to \infty\) (커널 안에 점이 충분히 많아야 평균이 적분에 수렴).

이 두 조건이 비모수 회귀의 표준 가정.

2.5 비유

큰 도시의 한 동네 (\(t\) 주위 작은 영역) 의 평균 인구 밀도를, 그 동네 안에 있는 모든 집의 거주자 수로 나눈 것. 동네가 충분히 작으면 (작은 \(h\)) 그 영역의 진짜 평균에 수렴, 동네 안 집이 충분히 많으면 (큰 \(Nh\)) 평균이 안정. \(\blacksquare\)


3 Problem 7.2: 식 (7.3) 의 정당화

3.1 문제

\(\mu\) 의 2 차 도함수가 유계라 가정하고 식 (7.3) \(F(h) = 2\mu(t) h + O(h^3)\) 를 정당화. 여기서 \(F(h) = \int_{t-h}^{t+h} \mu(x) \, dx\).

3.2 풀이

\(F(0) = 0\). \(F\) 의 미분:

\[ F'(h) = \mu(t + h) + \mu(t - h). \]

\(F'(0) = 2\mu(t)\).

\[ F''(h) = \mu'(t + h) - \mu'(t - h). \]

\(F''(0) = 0\) (대칭으로 자동).

\[ F'''(h) = \mu''(t + h) + \mu''(t - h). \]

\(F'''(0) = 2\mu''(t)\).

Taylor 전개 (h = 0 주위):

\[ F(h) = F(0) + F'(0) h + \frac{F''(0)}{2} h^2 + \frac{F'''(0)}{6} h^3 + O(h^5). \]

\(F(0) = F''(0) = 0\) 대입:

\[ F(h) = 2\mu(t) h + \frac{2 \mu''(t)}{6} h^3 + O(h^5) = 2\mu(t) h + \frac{\mu''(t)}{3} h^3 + O(h^5). \]

\(\mu''\) 가 유계 (\(|\mu''| \leq M\)) 면 \(|F(h) - 2\mu(t) h| \leq C h^3\). \(\blacksquare\)

3.3 직관: \(h^2\) 가 자동으로 사라지는 이유

박스 커널이 대칭 — \(K(x) = K(-x)\). 적분 \(F(h)\) 도 대칭 — \(t\) 양쪽에 같은 길이의 구간. 대칭 함수의 적분은 홀함수 (\(\mu'\) 의 차) 항이 자동 0.

이는 대칭 커널의 표준 결과 — 대칭 평활기는 항상 \(h^2\) 차 bias 를 자동 제거. 4 차 커널은 더 나아가 \(h^4\) 까지 제거.

3.4 비유

정사각형 영역의 평균값을 중심값으로 근사. 영역이 작으면 중심값과 평균값의 차이가 영역 크기의 거듭제곱으로 줄어듦. 정사각형의 대칭성이 1 차 항을 자동 0 으로 만들어 2 차 차이만 남김 (실제로는 3 차 — 대칭성 + 균일성).


4 Problem 7.3: 식 (7.5) 의 정당화

4.1 문제

식 (7.5) \(\text{Var}\left[\sum_m K\left(\frac{t - t_{nm}}{h}\right)(\varepsilon_n + \delta_{nm})\right] \approx M h (\tau^2 + \sigma^2) + M^2 h^2 \tau^2\) 를 정당화. \(h\) 에 어떤 조건이 필요한가?

4.2 풀이

박스 커널 가정. 단위 \(n\) 에 고정된 시점 \(\{t_{nm}\}\) 조건부 분산:

\[ \text{Var}\left[\sum_m K_h(t_{nm})(\varepsilon_n + \delta_{nm})\right] = \text{Var}\left[\sum_m K_h \varepsilon_n + \sum_m K_h \delta_{nm}\right], \]

\(K_h := K\left(\frac{t - t_{nm}}{h}\right)\).

두 항이 독립이므로 분산이 분리:

\[ = \underbrace{\text{Var}[\varepsilon_n \sum_m K_h]}_{(I)} + \underbrace{\sum_m K_h^2 \text{Var}[\delta_{nm}]}_{(II)}. \]

  1. \(K^2 = K\) (박스 커널), \(\text{Var}[\delta_{nm}] = \sigma^2\):

\[ (II) = \sigma^2 \sum_m K_h. \]

기댓값 \(E[\sum_m K_h] \approx M \cdot 2h \cdot 1 = 2Mh\) (시점이 균일이면 박스 안의 평균 점 수). 단순화로 \(\approx Mh\) (박스 폭 \(2h\), 단위 길이 \(1\) 가정).

  1. \(\varepsilon_n\) 가 시간 무관 상수이므로:

\[ (I) = \text{Var}[\varepsilon_n] \cdot \left(\sum_m K_h\right)^2 = \tau^2 \cdot (Mh)^2 = M^2 h^2 \tau^2. \]

  1. 에 추가로 \(\varepsilon_n\)\(\delta_{nm}\) 의 cross 항이 있다고 본문은 처리하지만, 본문 식에서는 단순하게 \(Mh(\tau^2 + \sigma^2) + M^2 h^2 \tau^2\) 로 표현 — 첫 항이 점별 분산의 합, 둘째 항이 같은 단위 내 cross 효과.

4.3 조건

  • \(h \to 0\).
  • \(Mh \to 0\) 또는 \(\to\) 상수 (박스 안의 점 수가 단위 내 다른 시점에 영향 안 미치도록).
  • 점근적으로 \(Mh \to 0\) 이지만 \(NMh \to \infty\) (전체 데이터 충분).

4.4 직관: 두 항의 의미 재확인

  • 첫 항 \(Mh(\tau^2 + \sigma^2)\) — 점별 분산의 합. 각 점의 독립 잡음 이 박스 안의 점 수만큼 누적.
  • 둘째 항 \(M^2 h^2 \tau^2\) — 같은 \(\varepsilon_n\) 가 박스 안의 모든 점에 영향. 단위 내 상관 이 곱셈적으로 누적.

이 둘째 항이 sparse FDA 분산 식의 핵심 — dense 비모수 회귀와의 결정적 차이.

4.5 비유: 측정 오차의 두 종류

100 개 측정의 평균을 낼 때:

  • 각 측정의 잡음 (첫 항) 은 \(\sqrt{100}\) 만큼 줄어들지만 (CLT),
  • 모든 측정에 공통인 시스템 오차 (둘째 항) 는 측정 횟수에 무관 — 시스템 자체를 바꿔야 줄어듦.

Sparse FDA 의 두 분산 항이 정확히 같은 패턴 — 측정 잡음 vs 단위 내 시스템 효과.


5 Problem 7.4: MSE 최적화

5.1 문제

Example 7.1.2 의 점근 MSE:

\[ \text{MSE}(h) = C h^2 + \frac{\tau^2 + \sigma^2}{NMh} + \frac{\tau^2}{N}. \]

  1. 미적분으로 \(h \geq 0\) 의 최적값 찾기. 최적 \(h\) 에서 bias 와 variance 의 크기 관계는?

  2. 최적 \(h\) 를 MSE 에 대입. \(M = N^\delta\) 일 때 MSE 의 점근 행동을 \(\delta\) 에 대해 분석.

(주: 본문 7.1 절은 Bias² ~ h^4 형태로 표시되었지만, Problem 7.4 는 Bias² 가 \(h^2\) 와 같은 형태인 MSE 식을 사용. 본문의 식을 직접 따른다.)

5.2 Part (a) 풀이

\(h\) 에 대한 미분:

\[ \frac{d \text{MSE}}{dh} = 2C h - \frac{\tau^2 + \sigma^2}{N M h^2}. \]

0 으로 놓고:

\[ 2 C h = \frac{\tau^2 + \sigma^2}{N M h^2} \implies h^3 = \frac{\tau^2 + \sigma^2}{2 C N M} \implies h^* = \left(\frac{\tau^2 + \sigma^2}{2 C N M}\right)^{1/3}. \]

따라서 \(h^* \sim (NM)^{-1/3}\).

(주: 본문 7.1 절의 결과 \((NM)^{-1/5}\) 는 Bias² ~ \(h^4\) 가정, 여기 Problem 7.4 는 Bias² ~ \(h^2\) 가정 — 두 식이 다른 환경.)

5.3 Bias 와 Variance 의 관계

최적 \(h^*\) 에서:

\[ \text{Bias}^2 = C h^{*2} = C \left(\frac{\tau^2 + \sigma^2}{2 C N M}\right)^{2/3}, \]

\[ \text{Var (h-항)} = \frac{\tau^2 + \sigma^2}{N M h^*} = (\tau^2 + \sigma^2) \left(\frac{2 C N M}{\tau^2 + \sigma^2}\right)^{1/3} \cdot \frac{1}{NM}. \]

대수적 정리:

\[ \frac{\text{Bias}^2}{\text{Var (h-항)}} = \frac{C h^{*2} \cdot NM h^*}{\tau^2 + \sigma^2} = \frac{C N M h^{*3}}{\tau^2 + \sigma^2}. \]

최적 \(h^{*3} = \frac{\tau^2 + \sigma^2}{2 C N M}\) 대입:

\[ \frac{\text{Bias}^2}{\text{Var (h-항)}} = \frac{C N M \cdot \frac{\tau^2 + \sigma^2}{2 C N M}}{\tau^2 + \sigma^2} = \frac{1}{2}. \]

따라서 최적 \(h\) 에서 Bias² = (1/2) × Variance(h-항). 즉 bias 가 variance 의 절반 — 표준 비모수 회귀 결과.

5.4 직관: Bias-Variance 균형의 정확한 비율

자주 “Bias² = Variance” 가 최적 균형이라 알려져 있으나, 정확한 비율은 손실 함수의 형태에 의존. MSE 손실에서는 Bias² : Variance = 1 : 2 (즉 Bias² = (1/2) Variance).

이는 미분의 구조에서 자동 도출 — Bias² ~ \(h^2\), Variance ~ \(h^{-1}\) 이면 \(\frac{d}{dh}\) 가 0 되는 점에서 자동으로 이 비율.

5.5 Part (b) 풀이

\(M = N^\delta\) 대입. \(h^* = (NM)^{-1/3} = (N \cdot N^\delta)^{-1/3} = N^{-(1+\delta)/3}\).

MSE 의 두 \(h\)-의존 항을 합치면 \(\sim h^2 \sim N^{-2(1+\delta)/3}\). \(h\)-무관 항은 \(N^{-1}\).

따라서:

\[ \text{MSE}^* \sim N^{-2(1+\delta)/3} + N^{-1}. \]

두 항 비교: \(-2(1+\delta)/3\) vs \(-1\).

\[ -\frac{2(1+\delta)}{3} \leq -1 \iff 2(1+\delta) \geq 3 \iff \delta \geq \frac{1}{2}. \]

5.6 시나리오

조건 MSE 행동
\(\delta > 1/2\) \(N^{-1}\) (모수적, \(h\)-무관 항이 지배)
\(\delta = 1/2\) 두 항 같은 차수
\(\delta < 1/2\) \(N^{-2(1+\delta)/3}\) (비모수적)

5.7 직관: 본문 7.1 과의 차이

본문 7.1 은 \((NM)^{-4/5}\)\(N^{-1}\) 을 비교하여 \(M \sim N^{1/4}\) 임계값. 여기서는 \(N^{-2(1+\delta)/3}\)\(N^{-1}\) 비교하여 \(\delta = 1/2\) 임계값 (즉 \(M \sim N^{1/2}\)).

차이는 Bias² 의 차수 가정 — 본문의 \(h^4\) vs Problem 7.4 의 \(h^2\). 더 매끄러운 가정 (\(\mu''\) 유계 + 대칭 커널 → \(h^4\)) 이 더 효율적인 추정 가능.

5.8 비유: 카메라 해상도와 셔터 속도

해상도 (대역폭의 정밀성) 와 셔터 속도 (분산 통제) 의 균형. 두 항 사이의 정확한 비율은 카메라의 광학 시스템 (모형의 매끄러움 가정) 에 의존 — Problem 7.4 는 한 가정, 본문 7.1 은 다른 가정에서의 균형.


6 Problem 7.5: 식 (7.9) 의 검증

6.1 문제

식 (7.9): \(\widehat{\boldsymbol{\beta}} = (\mathbf{Z}^T \mathbf{K}_h \mathbf{Z})^{-1} \mathbf{Z}^T \mathbf{K}_h \mathbf{Y}\) — Local polynomial regression 의 닫힌 해.

6.2 풀이

손실 함수:

\[ L(\boldsymbol{\beta}) = (\mathbf{Y} - \mathbf{Z}\boldsymbol{\beta})^T \mathbf{K}_h (\mathbf{Y} - \mathbf{Z}\boldsymbol{\beta}). \]

전개:

\[ L = \mathbf{Y}^T \mathbf{K}_h \mathbf{Y} - 2 \mathbf{Y}^T \mathbf{K}_h \mathbf{Z}\boldsymbol{\beta} + \boldsymbol{\beta}^T \mathbf{Z}^T \mathbf{K}_h \mathbf{Z} \boldsymbol{\beta}. \]

\(\boldsymbol{\beta}\) 에 대한 미분 (벡터 미적분):

\[ \frac{\partial L}{\partial \boldsymbol{\beta}} = -2 \mathbf{Z}^T \mathbf{K}_h \mathbf{Y} + 2 \mathbf{Z}^T \mathbf{K}_h \mathbf{Z} \boldsymbol{\beta}. \]

0 으로 놓고:

\[ \mathbf{Z}^T \mathbf{K}_h \mathbf{Z} \widehat{\boldsymbol{\beta}} = \mathbf{Z}^T \mathbf{K}_h \mathbf{Y} \implies \widehat{\boldsymbol{\beta}} = (\mathbf{Z}^T \mathbf{K}_h \mathbf{Z})^{-1} \mathbf{Z}^T \mathbf{K}_h \mathbf{Y}. \quad \blacksquare \]

6.3 직관: 가중 LS 의 표준 형태

이는 표준 가중 LS \(\widehat{\boldsymbol{\beta}} = (\mathbf{X}^T \mathbf{W} \mathbf{X})^{-1} \mathbf{X}^T \mathbf{W} \mathbf{Y}\) 의 특수 형태 — 가중 행렬 \(\mathbf{W} = \mathbf{K}_h\) (커널 가중).

이는 7.2 의 전체 framework — local polynomial = 각 시점에서의 가중 LS — 의 직접 확인.


7 Problem 7.6: 식 (7.10) 의 검증

7.1 문제

식 (7.10): \(\|f\|_{H_K}^2 = \sum_{j, k=1}^J \alpha_j \alpha_k K(s_j, s_k)\), where \(f \in A_K\) has the form \(f(t) = \sum_{j=1}^J \alpha_j K(t, s_j)\).

7.2 풀이

\(K\) 의 스펙트럼 분해 \(K(t, s) = \sum_i \lambda_i v_i(t) v_i(s)\), RKHS 내적의 정의:

\[ \langle f, g \rangle_{H_K} = \sum_i \frac{\int f v_i \cdot \int g v_i}{\lambda_i}. \]

\(f = \sum_j \alpha_j K(\cdot, s_j) = \sum_j \alpha_j \sum_i \lambda_i v_i(\cdot) v_i(s_j) = \sum_i \lambda_i \left(\sum_j \alpha_j v_i(s_j)\right) v_i(\cdot)\).

따라서 \(\int f v_i = \lambda_i \sum_j \alpha_j v_i(s_j)\).

\(f\) 자체 (\(f = g\)):

\[ \|f\|_{H_K}^2 = \sum_i \frac{\left(\lambda_i \sum_j \alpha_j v_i(s_j)\right)^2}{\lambda_i} = \sum_i \lambda_i \left(\sum_j \alpha_j v_i(s_j)\right)^2. \]

전개:

\[ = \sum_i \lambda_i \sum_{j, k} \alpha_j \alpha_k v_i(s_j) v_i(s_k) = \sum_{j, k} \alpha_j \alpha_k \sum_i \lambda_i v_i(s_j) v_i(s_k) = \sum_{j, k} \alpha_j \alpha_k K(s_j, s_k). \quad \blacksquare \]

7.3 직관: 핵으로 표현된 노름

이 식의 우아함:

RKHS 노름이 핵 자체의 값들의 이차 형식으로 표현 가능 — 스펙트럼 분해를 명시적으로 사용하지 않고도 계산.

이는 RKHS 의 실용적 장점 — 핵 함수만 정의하면 노름이 즉시 계산 가능.

7.4 비유: 거리의 두 표현

평면의 두 점 사이 거리 = \(\sqrt{(x_1 - x_2)^2 + (y_1 - y_2)^2}\) (좌표 기반) = \(\sqrt{r_1^2 + r_2^2 - 2 r_1 r_2 \cos\theta}\) (극좌표). 두 표현 모두 같은 거리, 단지 다른 좌표계.

RKHS 노름의 두 표현 — 스펙트럼 (\(\lambda_i, v_i\) 기반) vs 핵 직접 — 도 같은 패턴. 핵 직접 표현이 계산 편리.


8 Problem 7.7: RKHS 로 평균 추정

8.1 문제

Example 7.1.1 의 평균 함수를 RKHS Gaussian 핵 \(K(t, s) = \exp(-(t - s)^2)\) 로 추정. 벌점 모수 \(\lambda\) 를 CV 로 선택.

8.2 R 코드

library(stats)

# CATT 데이터의 (t, Y) 쌍이 있다고 가정
# t_obs, Y_obs: 관측 시점과 값의 벡터

# Gaussian 핵
gauss_K <- function(t, s) exp(-(t - s)^2)

# 핵 행렬
N <- length(t_obs)
K <- outer(t_obs, t_obs, gauss_K)

# 벌점 추정: alpha_hat = (K^T K + lambda K)^{-1} K^T Y
estimate_rkhs <- function(lambda) {
  alpha_hat <- solve(crossprod(K) + lambda * K, crossprod(K, Y_obs))
  list(alpha = alpha_hat, predict = function(t_new) {
    K_new <- outer(t_new, t_obs, gauss_K)
    K_new %*% alpha_hat
  })
}

# CV 로 lambda 선택 (subject 단위)
set.seed(2026)
n_folds <- 5
fold_id <- sample(1:n_folds, length(unique(subject_id)), replace = TRUE)
# (subject 단위 fold 할당)

lambda_grid <- 10^seq(-3, 3, length = 20)
cv_mse <- numeric(length(lambda_grid))

for (i in seq_along(lambda_grid)) {
  lambda <- lambda_grid[i]
  mse_folds <- sapply(1:n_folds, function(k) {
    train_idx <- which(fold_id[subject_id] != k)
    test_idx <- which(fold_id[subject_id] == k)

    K_train <- K[train_idx, train_idx]
    Y_train <- Y_obs[train_idx]

    alpha_hat <- solve(crossprod(K_train) + lambda * K_train,
                       crossprod(K_train, Y_train))

    K_test <- outer(t_obs[test_idx], t_obs[train_idx], gauss_K)
    Y_pred <- K_test %*% alpha_hat
    mean((Y_obs[test_idx] - Y_pred)^2)
  })
  cv_mse[i] <- mean(mse_folds)
}

best_lambda <- lambda_grid[which.min(cv_mse)]

# 최종 적합
final_fit <- estimate_rkhs(best_lambda)
t_grid <- seq(min(t_obs), max(t_obs), length = 200)
mu_hat <- final_fit$predict(t_grid)

plot(t_obs, Y_obs, pch = 19, cex = 0.5, col = rgb(0, 0, 0, 0.3),
     xlab = "Week", ylab = "VAS",
     main = paste("RKHS Gaussian (lambda =", round(best_lambda, 4), ")"))
lines(t_grid, mu_hat, col = "red", lwd = 2)

8.3 직관: CV 의 subject-level 분할

7.2 에서 강조한 원칙 — sparse FDA 의 CV 는 subject 단위 분할. 한 환자의 점들이 같은 fold 에 속해야 단위 내 상관이 fold 평가에 영향 없음.

이 원칙을 무시하면 CV 가 under-smoothing — 너무 작은 \(\lambda\) 선택.


9 Problem 7.8: Baseline 빼고 gam 적합

9.1 문제

Example 7.1.1 의 평균 함수를 gam 으로 적합하되 baseline 을 제외하고 다른 관측에서 baseline 을 차감. 이 처리가 CV 의 합리적 모수 선택에 도움이 되는가?

9.2 풀이 스케치

CATT 데이터의 baseline (treatment 시작 전 측정) 과 후속 측정 사이에 큰 jump 가 있다 — 치료 시작 직후 빠른 회복. 이 불연속이 매끄러움 가정과 모순.

대안: 각 환자의 baseline 을 빼서 정규화:

Y_centered <- Y_obs - rep(baseline_per_subject, times = M_per_subject)

이러면 모든 환자가 baseline 0 에서 시작 — 불연속 제거. 결과 CV 가 더 합리적 \(\lambda\) 선택.

9.3 직관: 데이터의 구조에 맞는 전처리

자동 알고리즘이 문제를 잘 해결하지 못하면, 데이터의 구조 (불연속) 를 모형이 더 잘 처리할 수 있도록 전처리. CATT 의 baseline 정규화가 그 예시.

이는 통계 분석의 일반 원칙 — 모형의 가정에 맞도록 데이터 변환.


10 Problem 7.9: 식 (7.7) 이 (7.6) 을 최소화

10.1 문제

식 (7.6) \(\widehat{\mu}(t) = \arg\min_\alpha \sum_n \sum_m K_h (Y_{nm} - \alpha)^2\) 의 해가 식 (7.7) 의 NW 추정량임을 보임.

10.2 풀이

손실의 \(\alpha\) 에 대한 미분:

\[ \frac{d}{d\alpha} \sum_n \sum_m K_h (Y_{nm} - \alpha)^2 = -2 \sum_n \sum_m K_h (Y_{nm} - \alpha). \]

0 으로:

\[ \sum K_h Y_{nm} = \alpha \sum K_h \implies \widehat{\alpha} = \frac{\sum K_h Y_{nm}}{\sum K_h}. \]

이는 식 (7.7) 의 NW 추정량. \(\blacksquare\)

10.3 직관: 가중 평균이 가중 LS 의 해

NW 추정량은 가중 평균 — 각 점에 커널 가중치. 이는 가중 LS 의 가장 단순한 형태 (\(P = 0\) local polynomial).


11 Problem 7.10: CATT 공분산 추정

11.1 R 풀이 (코드 스케치)

library(mgcv)

# 잔차 형성
mu_hat <- predict(mu_fit, newdata = data.frame(t = t_obs))
residuals <- Y_obs - mu_hat

# 비대각 cross product 만 사용
cross_data <- data.frame()
for (n in unique(subject_id)) {
  idx <- which(subject_id == n)
  if (length(idx) < 2) next
  pairs <- combn(length(idx), 2)
  for (k in 1:ncol(pairs)) {
    i1 <- idx[pairs[1, k]]
    i2 <- idx[pairs[2, k]]
    cross_data <- rbind(cross_data,
      data.frame(t1 = t_obs[i1], t2 = t_obs[i2],
                 prod = residuals[i1] * residuals[i2]))
    cross_data <- rbind(cross_data,   # 대칭 추가
      data.frame(t1 = t_obs[i2], t2 = t_obs[i1],
                 prod = residuals[i1] * residuals[i2]))
  }
}

# Tensor product spline 으로 이변량 평활
cov_fit <- gam(prod ~ te(t1, t2, k = c(8, 8)), data = cross_data, method = "REML")

# 격자 위에서 평가
t_grid <- seq(0, 100, length = 50)
grid <- expand.grid(t1 = t_grid, t2 = t_grid)
c_hat <- matrix(predict(cov_fit, newdata = grid), 50, 50)

# 양정치 보정
eig <- eigen(c_hat, symmetric = TRUE)
eig$values[eig$values < 0] <- 0
c_hat_psd <- eig$vectors %*% diag(eig$values) %*% t(eig$vectors)

# 잡음 분산
diag_smooth <- gam(I(residuals^2) ~ s(t_obs, k = 20), method = "REML")
c_tilde_diag <- predict(diag_smooth, newdata = data.frame(t_obs = t_grid))
sigma2_hat <- pmax(c_tilde_diag - diag(c_hat_psd), 0)

11.2 직관

이는 7-2 포스트의 R 구현과 동일. CATT 데이터의 표준 sparse FDA 분석 워크플로우.


12 Problem 7.11: Representer Theorem 증명 (x ∈ A_K 경우)

12.1 문제

\(x \in A_K\) (즉 \(x(t) = \sum_j \alpha_j K(t, s_j)\) 형태) 에 대해 Theorem 7.2.1 (\(x(t) = \langle x, K_t \rangle_{H_K}\)) 증명.

12.2 풀이

\(x = \sum_j \alpha_j K_{s_j}\) 라 표기 (\(K_{s_j}(s) = K(s_j, s) = K(s, s_j)\)).

내적:

\[ \langle x, K_t \rangle_{H_K} = \langle \sum_j \alpha_j K_{s_j}, K_t \rangle = \sum_j \alpha_j \langle K_{s_j}, K_t \rangle. \]

\(\langle K_{s_j}, K_t \rangle\) 의 계산. RKHS 내적의 정의:

\[ \langle K_{s_j}, K_t \rangle = \sum_i \frac{\int K_{s_j}(u) v_i(u) \, du \cdot \int K_t(u) v_i(u) \, du}{\lambda_i}. \]

\(K_t(u) = \sum_i \lambda_i v_i(t) v_i(u)\) 로 (스펙트럼 분해), \(\int K_t(u) v_i(u) \, du = \lambda_i v_i(t)\) (정규직교성).

같은 방식으로 \(\int K_{s_j}(u) v_i(u) \, du = \lambda_i v_i(s_j)\).

대입:

\[ \langle K_{s_j}, K_t \rangle = \sum_i \frac{\lambda_i v_i(s_j) \cdot \lambda_i v_i(t)}{\lambda_i} = \sum_i \lambda_i v_i(s_j) v_i(t) = K(s_j, t). \]

따라서:

\[ \langle x, K_t \rangle = \sum_j \alpha_j K(s_j, t) = \sum_j \alpha_j K(t, s_j) = x(t). \quad \blacksquare \]

12.3 직관: 핵의 두 가지 역할

\(K(t, s)\) 가:

  1. 함수 표현의 building block\(x(t) = \sum_j \alpha_j K(t, s_j)\).
  2. 점 평가의 표현 함수\(K_t\) 가 점 평가 범함수의 RKHS 표현.

두 역할이 일치 (\(K_t\)\(K\) 의 한 슬라이스) — 이것이 RKHS 의 우아함.

12.4 비유: 같은 도구의 두 사용

해머가 (1) 못을 박는 도구 + (2) 못을 빼는 도구의 두 역할. 같은 도구지만 사용 방식이 다름. 핵 \(K\) 도 (1) 함수 구성 + (2) 점 평가의 두 역할 — 같은 핵의 두 사용.


13 Problem 7.12: RKHS 최소화의 형태

13.1 문제

벌점 RKHS 최소화의 해가 \(\widehat{\mu}(t) = \sum_{n, m} \widehat{\alpha}_{nm} K(t, t_{nm})\) 형태임을 증명.

13.2 풀이 (Representer Theorem 의 응용)

\(H_K\)\(V \oplus V^\perp\) 로 분해. \(V = \text{span}\{K_{t_{nm}}: n, m\}\), \(V^\perp\) = 직교 보충.

임의의 \(\mu \in H_K\)\(\mu = \mu_V + \mu_{V^\perp}\) 로 분해.

관측에서의 평가:

\[ \mu(t_{nm}) = \langle \mu, K_{t_{nm}} \rangle = \langle \mu_V, K_{t_{nm}} \rangle + \underbrace{\langle \mu_{V^\perp}, K_{t_{nm}} \rangle}_{= 0}. \]

\(K_{t_{nm}} \in V\) 이고 \(\mu_{V^\perp} \in V^\perp\) 이므로 두 번째 항은 0.

따라서 \(\mu(t_{nm}) = \mu_V(t_{nm})\)잔차 제곱합은 \(V\) 부분에만 의존.

노름:

\[ \|\mu\|^2 = \|\mu_V\|^2 + \|\mu_{V^\perp}\|^2 \geq \|\mu_V\|^2. \]

따라서 \(\mu_{V^\perp} \neq 0\) 이면 잔차 제곱합은 같지만 노름은 증가 — 손실 함수가 더 큼. 최소화하려면 \(\mu_{V^\perp} = 0\), 즉 \(\mu \in V\).

따라서 최소화 해의 형태:

\[ \widehat{\mu}(t) = \sum_{n, m} \widehat{\alpha}_{nm} K(t, t_{nm}). \quad \blacksquare \]

13.3 직관: 벌점이 단순화하는 이유

벌점 (RKHS 노름) 이 없으면 임의의 \(\mu_{V^\perp}\) 추가가 가능 — 무한차원 자유도. 벌점이 \(\mu_{V^\perp} = 0\) 을 강제 — 추정해야 할 차원이 데이터 점 수 (\(\sum M_n\)) 로 환원.

이는 RKHS 의 우아한 결과 — 무한차원 최적화가 유한차원으로 자동 환원.

13.4 비유: 가장 짧은 경로

A 에서 B 로 가는 가장 짧은 경로는 직선. 우회로를 추가하면 (\(V^\perp\) 방향) 더 멀어짐. 벌점이 “거리 최소화” 역할이므로 자동으로 직선 경로 (\(V\) 방향만) 선택.


14 Problem 7.13: Canadian Weather 의 RKHS 평균 (두 핵)

14.1 문제

Canadian Weather 의 평균 함수를 RKHS 로 추정. 두 핵 비교:

  1. Gaussian: \(K(t, s) = \exp(-(t-s)^2)\)
  2. Periodic: \(K(t, s) = \exp(-\sin^2(\pi |t - s|))\)

14.2 R 풀이 스케치

library(fda)
data(CanadianWeather)

# 일평균 기온 (35 관측소 × 365 일)
temps <- CanadianWeather$dailyAv[, , "Temperature.C"]
day <- 1:365 / 365   # 도메인 [0, 1] 로 재척도

# Long format
long_data <- data.frame(
  station = rep(1:35, each = 365),
  t = rep(day, 35),
  Y = as.vector(temps)
)

# (a) Gaussian RKHS
gauss_K <- function(t, s) exp(-(t - s)^2)
K_gauss <- outer(long_data$t, long_data$t, gauss_K)
alpha_gauss <- solve(crossprod(K_gauss) + lambda * K_gauss,
                     crossprod(K_gauss, long_data$Y))

# (b) Periodic RKHS
periodic_K <- function(t, s) exp(-sin(pi * abs(t - s))^2)
K_periodic <- outer(long_data$t, long_data$t, periodic_K)
alpha_periodic <- solve(crossprod(K_periodic) + lambda * K_periodic,
                        crossprod(K_periodic, long_data$Y))

14.3 결과 비교 직관

  • Gaussian: 일반 매끄러움. 1 월 (저점) → 7 월 (정점) → 12 월 (저점) 의 부드러운 곡선.
  • Periodic: 주기 구조 강제. 1 월 1 일과 12 월 31 일이 자동으로 비슷한 값 (도메인의 양 끝이 가까운 것으로 처리).

기온 데이터는 본질적으로 주기적 — Periodic 핵이 더 자연스러운 결과.

14.4 직관: 도메인 지식이 핵 선택을 결정

같은 데이터에 다른 핵을 적용하면 다른 결과. 도메인 지식 (주기성) 을 반영하는 핵이 더 적절.

이것이 RKHS 의 강점 — Spline 이나 local polynomial 같은 일반 도구는 주기성을 직접 강제하기 어렵지만, 핵 선택으로 자동 처리.

14.5 비유: 적합한 옷

같은 사람이라도 상황 (직장·운동·결혼식) 에 맞는 옷이 있음. 데이터의 본성 (주기·매끄러움·경계 조건) 에 맞는 핵이 다른 결과를 낸다.


15 Problem 7.14: Canadian Weather 공분산의 양정치성

15.1 문제

Problem 7.13 의 데이터로 공분산 함수 추정. 격자에서 평가하여 대칭성, 양정치성 확인.

15.2 풀이 스케치

평활 추정 \(\widehat{c}(t, s)\) 가:

  • 대칭 — 평활이 \((t, s)\)\((s, t)\) 를 동일하게 처리하면 자동 대칭. 비대칭이면 명시적 대칭화: \(\widehat{c}^{sym}(t, s) = (\widehat{c}(t, s) + \widehat{c}(s, t)) / 2\).
  • 양정치 보정 — 7-2 의 표준 절차. 음의 고유값을 0 으로.

15.3 직관

7-2 의 양정치 보정 절차의 구체적 적용. Canadian Weather 같은 매끄러운 데이터에서는 음의 고유값이 작거나 없을 가능성 — 공분산 평활이 자연스럽게 양정치.


16 Problem 7.15: BLUP of Y(t)

16.1 문제

\(\{Y(t): t \in \mathcal{T}\}\) 가 가우스 과정 (평균 \(\mu(t)\), 공분산 \(c(t, s) + \sigma^2(t) \delta_{t s}\)). \(Y(t_1), \ldots, Y(t_m)\) 관측. \(Y(t)\) 의 BLUP 은?

16.2 풀이

PACE 의 일반화. 결합 가우스:

\[ \begin{pmatrix} Y(t) \\ \mathbf{Y} \end{pmatrix} \sim N\left(\begin{pmatrix} \mu(t) \\ \boldsymbol{\mu} \end{pmatrix}, \boldsymbol{\Sigma}\right), \]

\(\boldsymbol{\Sigma} = \begin{pmatrix} c(t, t) + \sigma^2(t) & \mathbf{c}_{12}^T \\ \mathbf{c}_{12} & \boldsymbol{\Sigma}_{22} \end{pmatrix}\), \(\mathbf{c}_{12} = (c(t, t_1), \ldots, c(t, t_m))^T\), \(\boldsymbol{\Sigma}_{22} = [c(t_i, t_j) + \sigma^2(t_i) \delta_{ij}]\).

조건부 기댓값 (BLUP):

\[ \widehat{Y}(t) = \mu(t) + \mathbf{c}_{12}^T \boldsymbol{\Sigma}_{22}^{-1} (\mathbf{Y} - \boldsymbol{\mu}). \quad \blacksquare \]

16.3 직관: PACE 의 곡선 직접 BLUP

7-2 의 PACE 가 점수 BLUP. 여기는 곡선 자체의 BLUP — 같은 framework, 다른 출력.

가우스 가정 하 두 BLUP 은 동등 (\(Y(t) = \sum_j \xi_j v_j(t)\) + 잡음, 점수의 BLUP 으로 \(Y(t)\) 재구성). 그러나 곡선 BLUP 이 더 직접 — FPCA 차원 축소 없이 바로 곡선 예측.

실무 차이: 점수 BLUP 은 차원 축소 + 후속 분석에 유용, 곡선 BLUP 은 평활 곡선 시각화에 직접.

16.4 비유: 사진의 두 복원 방식

저해상도 사진의 고해상도 복원 두 방식:

  1. 특징 추출 후 재구성 (PCA 기반) — 얼굴 특징을 추출 후 합성.
  2. 픽셀 직접 복원 (super-resolution) — 픽셀 값을 직접 예측.

두 방식이 본질적으로 같은 정보를 사용하지만 표현 방식이 다름. PACE 점수 BLUP vs 곡선 BLUP 도 같은 패턴.


17 Problem 7.16: Local Polynomial 과 도함수

17.1 문제

식 (7.8) 의 local polynomial 회귀와 Taylor 전개 비교. (a) \(\mu^{(i)}(t)\)\(\beta_i\) 의 관계는? (b) 시뮬레이션으로 검증 (\(\mu(t) = \sin(2\pi t)\), \(P = 1\), \(h = N^{-1/5}\)).

17.2 Part (a) 풀이

Taylor 전개:

\[ \mu(t) = \sum_{i=0}^P \frac{\mu^{(i)}(t_0) (t - t_0)^i}{i!} + o(|t - t_0|^P). \]

Local polynomial 모형 (7.8):

\[ \mu(t) \approx \sum_{i=0}^P \beta_i (t - t_0)^i. \]

두 식 비교:

\[ \boxed{ \beta_i = \frac{\mu^{(i)}(t_0)}{i!}. } \]

17.3 직관: 도함수 추정의 자동 부산물

Local polynomial 회귀의 계수 \(\widehat{\beta}_i\) 가 자동으로 도함수의 추정. 즉:

  • \(\widehat{\beta}_0 = \widehat{\mu}(t_0)\) — 함수 값.
  • \(\widehat{\beta}_1 = \widehat{\mu}'(t_0)\) — 1 차 도함수.
  • \(\widehat{\beta}_2 = \widehat{\mu}''(t_0) / 2\) — 2 차 도함수의 절반 (계수에 \(1/2!\)).
  • \(\widehat{\beta}_i = \widehat{\mu}^{(i)}(t_0) / i!\)\(i\) 차 도함수.

17.4 Part (b) 시뮬레이션

set.seed(2026)
N <- 1000; M <- 5
mu_f <- function(t) sin(2 * pi * t)

# 데이터 생성
t_all <- numeric(0); Y_all <- numeric(0)
for (n in 1:N) {
  t_n <- runif(M)
  eps_n <- rnorm(1, sd = 1)        # tau = 1
  delta_n <- rnorm(M, sd = sqrt(0.5))  # sigma^2 = 0.5
  Y_n <- mu_f(t_n) + eps_n + delta_n
  t_all <- c(t_all, t_n); Y_all <- c(Y_all, Y_n)
}

# Local linear (P = 1) with h = N^{-1/5}
h <- (N * M)^(-1/5)

# 50 격자 점에서 평가
t_grid <- seq(0, 1, length = 50)
mu_hat <- numeric(50)
mu_prime_hat <- numeric(50)

for (i in seq_along(t_grid)) {
  t0 <- t_grid[i]
  weights <- ifelse(abs(t_all - t0) < h, 1, 0)
  if (sum(weights) < 3) next
  fit <- lm(Y_all ~ I(t_all - t0), weights = weights)
  mu_hat[i] <- coef(fit)[1]
  mu_prime_hat[i] <- coef(fit)[2]   # = mu'(t0) (since 1!=1)
}

# 진짜 함수와 비교
par(mfrow = c(1, 2))
plot(t_grid, mu_hat, type = "l", col = "red", lwd = 2,
     xlab = "t", ylab = expression(mu(t)),
     main = "mu(t) estimate")
lines(t_grid, mu_f(t_grid), col = "black", lty = 2, lwd = 2)
legend("topright", c("Estimate", "True"), col = c("red", "black"), lty = c(1, 2))

plot(t_grid, mu_prime_hat, type = "l", col = "red", lwd = 2,
     xlab = "t", ylab = expression(mu*minute*(t)),
     main = "mu'(t) estimate")
lines(t_grid, 2*pi*cos(2*pi*t_grid), col = "black", lty = 2, lwd = 2)
legend("topright", c("Estimate", "True"), col = c("red", "black"), lty = c(1, 2))

17.5 직관: Local polynomial 의 두 가지 출력

이 결과의 매력: 하나의 적합으로 함수와 도함수 동시 추정. Local linear 적합 시 절편 = 함수, 기울기 = 1 차 도함수. Local quadratic 적합 시 추가로 2 차 도함수.

이는 도함수 추정의 다른 도구 (수치 미분 + 평활) 보다 더 효율적이고 안정적.

17.6 비유: 파일럿의 다중 계기판

비행기의 한 GPS 모듈이 위치 + 속도 + 가속도를 동시에 출력. 위치만 측정한 후 차분으로 속도 계산하면 잡음 증폭. 동시 추정이 더 정확.

Local polynomial 의 도함수 추정도 같은 원리 — 별도 차분보다 더 안정적.


18 Problem 7.17: Gaussian vs Exponential 핵의 매끄러움

18.1 문제

  1. Gaussian 핵 RKHS 의 함수가 적어도 1 회 미분 가능함을 보임.
  2. Exponential 핵 RKHS 의 함수가 미분 가능 보장이 안 되는 이유.

18.2 Part (a) 풀이 스케치

Real analysis 의 표준 결과 사용:

Theorem: 함수 수열 \(\{f_n\}\) 이 점별 수렴 + 도함수 \(\{f_n'\}\) 이 균등 수렴이면, 극한 함수 \(f\) 의 도함수 = 도함수의 극한.

Gaussian 핵 RKHS 의 임의 함수:

\[ x(t) = \sum_j \alpha_j K(t, s_j) = \sum_j \alpha_j e^{-\sigma (t - s_j)^2}. \]

각 항 \(e^{-\sigma(t-s_j)^2}\) 는 무한 미분 가능 + 도함수 \(-2\sigma(t-s_j) e^{-\sigma(t-s_j)^2}\) 도 매끄럽다.

부분합 \(x_J = \sum_{j \leq J} \alpha_j K(t, s_j)\) 와 그 도함수가 균등 수렴 (RKHS 노름 유한 조건 하). Theorem 에 의해 \(x\) 가 미분 가능.

18.3 Part (b) 풀이

Exponential 핵:

\[ K(s, s') = e^{-\sigma|s - s'|}. \]

\(|s - s'|\) 의 도함수가 \(s = s'\) 에서 정의되지 않음 (절댓값의 도함수 = 부호 함수, \(s = s'\) 에서 점프).

따라서 \(K(\cdot, s')\)\(s = s'\) 에서 미분 가능하지 않음. 그 핵의 선형 결합으로 만들어진 함수도 일반적으로 미분 가능 보장 안 됨.

18.4 직관: 핵의 매끄러움이 RKHS 의 매끄러움

\(K(s, s')\)\(s = s'\) 에서 얼마나 매끈한지가 RKHS 의 함수 매끄러움 결정.

\(s = s'\) 의 매끄러움 RKHS 의 함수
Gaussian (\(e^{-\sigma(s-s')^2}\)) 무한 미분 가능 무한 미분 가능
Sobolev (m차) \(m\) 회 미분 가능 \(m\) 회 미분 가능
Exponential (\(e^{-\sigma|s-s'|}\)) 연속이지만 1 차 미분 점프 연속이지만 미분 안 될 수 있음

이는 Karhunen-Loève 이론의 일반 결과 — 가우스 과정의 sample path 매끄러움이 공분산의 대각 매끄러움에 의해 결정.

18.5 비유: 줄의 탄성

핵이 함수 공간의 “탄성” 을 결정.

  • Gaussian = 매우 탄력적인 줄. 작은 흔들림도 부드러운 곡선.
  • Exponential = 덜 탄력적. 흔들림이 거친 모양 만듦.

도구 (핵) 의 본성이 결과 (RKHS 의 함수) 의 본성을 결정.


19 17 문제의 통합 시각

19.1 한 줄 요약

Ch.7 의 17 개 문제는 sparse FDA 의 점근 분석 (7.1~7.4: bias-variance 분해, MSE 최적화, 임계값 시나리오), local polynomial 과 RKHS 의 닫힌 해 (7.5~7.6, 7.11~7.12), R 코드 응용 (7.7~7.10, 7.13), 양정치 보정·BLUP·도함수 추정 (7.14~7.16), 핵 매끄러움 (7.17) 을 통합 검증한다. 특히 Problem 7.4 의 MSE 최적화와 Problem 7.16 의 Taylor 전개 도함수 추정이 sparse FDA 의 핵심 통찰이다.

19.2 그룹별 정리

그룹 문제 핵심 도구
점근 분석 7.1~7.4 Bias 식 (7.2), Taylor 전개 (7.3), 분산 (7.5), MSE 최적화
닫힌 해 7.5~7.6, 7.11~7.12 Local poly (7.9), RKHS 노름 (7.10), Representer Theorem
R 응용 7.7~7.10, 7.13 RKHS Gaussian, gam, CATT 공분산, Canadian Weather
보정·BLUP 7.14~7.16 양정치 보정, 곡선 BLUP, 도함수 동시 추정
핵 매끄러움 7.17 Gaussian (∞ 미분) vs Exponential (미분 불보장)

19.3 본문과의 매핑

본문 절 관련 문제
7.1 (도입 + 점근) 7.1, 7.2, 7.3, 7.4 (점근 분석 모두)
7.2 (평균 추정) 7.5 (local poly), 7.6 (RKHS 노름), 7.7~7.9 (R), 7.11~7.13 (RKHS), 7.16, 7.17
7.3 (공분산) 7.10 (CATT), 7.14 (양정치)
7.4 (PACE) 7.15 (곡선 BLUP)
7.5 (sparse 회귀) (없음 — 7.5 의 식은 본문 자체에서 충분)

연습문제 7.4 의 MSE 최적화가 7.1 의 정확한 검증, 7.16 의 도함수 추정이 7.2 local polynomial 의 활용 확장. 17 문제가 Ch.7 본문의 모든 핵심을 검증·확장.


20 관련 주제

선행 지식

후속 주제

관련 개념

Subscribe

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