제1부 수학적 기초

1.21 확률·통계 III: 확률 과정·가우시안 과정·몬테카를로 방법

1.21.1 이 장에서 배우는 것 / 선수 지식 / 특히 중요한 독자

이 장은 시간 또는 입력에 따라 변하는 불확실성을 다루는 세 가지 도구, 곧 확률 과정, 가우시안 과정, 몬테카를로 방법을 정리한다. 1.19은 한 시점의 확률변수를, 1.20은 모수의 추정을 다룬다. 이 장은 확률변수의 열과 함수에 대한 분포, 그리고 닫힌 형식으로 계산할 수 없는 기댓값을 표본으로 계산하는 방법을 다룬다. 센서 잡음의 시간 상관, 필터의 프로세스 잡음(9.3), 파티클 필터(9.5), 표본 기반 제어(12.19)가 이 도구에 의존한다. 다룰 질문은 다섯 가지다. 정상성·자기상관·전력 스펙트럼 밀도는 어떻게 이어지며 위너 과정과 오른스타인-울렌벡(OU) 과정은 무엇인가. 가우시안 과정은 커널로 무엇을 정하며 사후 예측·주변 우도·비용은 어떠한가. 표본 기반 기댓값의 오차는 어떻게 줄고 기각·중요도 표본추출의 효율은 어떤 조건에서 급감하는가. 메트로폴리스-헤이스팅스의 제안 폭은 효율을 어떻게 정하는가. R-hat과 유효 표본 수는 무엇을 검사하지 못하는가.

배우는 것:

• 확률 과정 — 정상성, 자기공분산, 전력 스펙트럼 밀도, 위너 과정, 화이트 노이즈, 마르코프성, OU 과정(AR(1)의 정확한 이산화).
• 가우시안 과정 — 평균·공분산 함수, RBF·Matérn 커널, 사후 예측, 주변 우도, 하이퍼파라미터, $O(N^3)$ 비용, 희소 근사(SR), 칼만 필터와의 등가.
• 몬테카를로 — $N^{-1/2}$ 오차, 기각 표본추출, 중요도 표본추출, 유효 표본 수와 가중치 퇴화, 분산 감소.
• 마르코프 체인 몬테카를로 — 메트로폴리스-헤이스팅스, 수락률, 적분 자기상관 시간, 분할 R-hat과 그 한계.
• 라이브러리 — scikit-learn, SciPy, NumPy, ArviZ(arviz-stats)와의 비교.

선수 지식: 1.19(가우시안·조건부·베이즈 갱신), 1.2(촐레스키·최소제곱), 1.1(LU), 1.25(푸리에 변환·주파수 응답). 추정 이론은 1.20, 구간 해석은 1.22, 필터는 9.1·9.3·9.5, 베이지안 최적화·MPPI는 1.16가 다룬다.

특히 중요한 독자: 센서 잡음 모델과 필터의 잡음 공분산을 정하는 엔지니어, 가우시안 과정 기반 동역학·궤적 모델을 구현하는 연구자, 표본 기반 계획·제어·평가 결과의 통계적 신뢰도를 판단하는 개발자.

출처 방침. 서지는 Crossref 응답으로 확인했다(목록은 이 장 끝). 직접 읽은 자료: Rasmussen·Williams 『Gaussian Processes for Machine Learning』의 공개 PDF(2.2·알고리즘 2.1·2.3·4.2·5.4.1·8장 도입·8.3.1·8.3.3·부록 B.2.1), Vehtari 외(2021)의 §1·§3.1·§3.2 일부와 권고 부분, Martino 외의 초록·§1~3, Owen의 몬테카를로 교재 9.2~9.3, JCGM 101:2008의 §1·§7.2, scikit-learn 1.9.1 소스(_gpr.py, kernels.py의 WhiteKernel)와 설치된 arviz-stats 1.3.3의 base/diagnostics.py. 직접 실행한 것: scikit-learn 1.9.1, SciPy 1.18.1, NumPy 2.5.4, arviz-stats 1.3.3(Python 3.12.3). 읽지 못한 것: Metropolis 외(1953)·Hastings(1970)·Gelman·Rubin(1992)·Geyer(1992)·Roberts 외(1997)·Doob(1942)·Wiener(1930)·Khintchine(1934) 등 원 논문의 본문(서지만 확인), Hartikainen·Särkkä(2010), Särkkä·Solin의 교재, Titsias와 Quiñonero-Candela·Rasmussen의 본문. 이 장의 수치는 모두 순수 파이썬 구현으로 계산했으며 1.2의 촐레스키와 1.19의 가우시안 표본 추출을 불러와 사용했다.

1.21.2 개념 — 시간 상관·함수의 분포·표본 계산

확률변수가 한 시점의 값을 다루는 대상이라면 확률 과정은 시간으로 색인된 확률변수의 모음이다. 센서 잡음이 시간에 따라 상관되면 표본이 독립이라는 가정이 깨지고, 평균의 정확도는 표본 수가 아니라 상관 시간으로 정해지는 유효 표본 수에 따른다. 가우시안 과정은 색인을 시간에서 임의의 입력(상태·위치)으로 넓혀 함수에 분포를 둔다. 유한 개 입력에서의 값이 결합 가우시안이므로 1.19의 조건부 분포 공식이 그대로 회귀가 된다. 몬테카를로 방법은 기댓값을 표본 평균으로 대신하며, 표본이 독립이면 오차는 $N^{-1/2}$로 줄고 표본이 상관되면(마르코프 체인) 유효 표본 수로 정해진다.

세 도구는 서로 이어진다. OU 과정은 가우시안 과정의 특수한 경우(Matérn $\nu=1/2$ 커널)이고, 이 커널의 시간 입력 회귀는 칼만 필터와 같은 계산이다. 자기상관으로 정의한 유효 표본 수 $N/(1+2\sum\rho_t)$는 마르코프 체인의 효율과 상관 센서 잡음의 평균 정확도에 같은 식으로 적용된다. 아래 도해는 OU 과정의 자기상관이다.

1.21.3 원리

1.21.3.1 확률 과정, 정상성, 전력 스펙트럼 밀도

결론: 광의 정상 과정의 자기공분산과 전력 스펙트럼 밀도는 푸리에 변환 쌍이다. 상관된 표본으로 평균을 추정하면 오차는 $N^{-1/2}$가 아니라 유효 표본 수 $N_{\mathrm{eff}}$에 따라 줄어든다.

확률 과정 $\{x(t)\}$의 평균은 $m(t)=\mathbb E[x(t)]$, 자기공분산은 $R(t,s)=\mathrm{Cov}[x(t),x(s)]$이다. 평균이 상수이고 $R(t,s)=R(t-s)$이면 광의 정상(wide-sense stationary)이다. 모든 유한 차원 분포가 시간 이동에 불변이면 엄밀 정상이며, 가우시안 과정은 평균과 공분산이 분포를 정하므로 두 정의가 일치한다. 한 경로의 시간 평균이 앙상블 평균에 수렴하는 에르고드성은 단일 기록으로 통계량을 추정하는 근거가 되는 가정이다. 자기상관은 $\rho(\tau)=R(\tau)/R(0)$이고 전력 스펙트럼 밀도(PSD)는 자기공분산의 푸리에 변환이다(위너-힌친 정리, Wiener 1930; Khintchine 1934; 정리 번호는 (확인 필요)).

$$S(\omega)=\int_{-\infty}^{\infty}R(\tau)\,e^{-i\omega\tau}\,d\tau,\qquad R(0)=\frac{1}{2\pi}\int_{-\infty}^{\infty}S(\omega)\,d\omega$$

이 장은 $\omega$[rad/s]에 대한 양쪽 밀도를 사용한다. 데이터시트의 잡음 밀도($\sqrt{\mathrm{Hz}}$당 신호 단위)는 한쪽 밀도의 제곱근이며, 주파수 $f$[Hz]에 대한 한쪽 밀도는 $2S(2\pi f)$이다(1.21.5⑤에서 scipy.signal.welch와 확인). 자기상관이 $\rho_k$인 $N$개 표본의 평균은 $\mathrm{Var}(\bar x)\approx\frac{R(0)}{N}\bigl(1+2\sum_{k\ge1}\rho_k\bigr)$이므로 유효 표본 수는 $N_{\mathrm{eff}}=N/(1+2\sum_{k\ge1}\rho_k)$이다. 이 식은 마르코프 체인의 유효 표본 수와 같다(Vehtari 외 2021, 식 7). 화이트 노이즈는 $R(\tau)=q\,\delta(\tau)$, $S(\omega)=q$로 평탄하며 분산이 무한하므로 물리적 신호가 아니다. 간격 $\Delta t$로 평균한 표본은 분산이 $q/\Delta t$인 이산 화이트 잡음에 대응한다.

1.21.3.2 위너 과정, 마르코프성, OU 과정

결론: 위너 과정은 독립 증분을 가진 비정상 가우시안 과정이고 화이트 노이즈의 적분이다. 정상 가우시안-마르코프 과정은 지수 자기상관을 가진 OU 과정이며, 일정 간격 표본은 AR(1)이다.

위너 과정은 $W(0)=0$이고 증분 $W(t)-W(s)\sim\mathcal N(0,t-s)$가 서로 독립이며 $\mathrm{Cov}[W(t),W(s)]=\min(t,s)$이다. 분산이 $t$에 비례해 커지므로 비정상이고, 경로는 연속이나 어디에서도 미분되지 않는 것으로 알려져 있다 (확인 필요: 정리 번호). 센서 바이어스의 랜덤 워크는 이 모델로 표현되며 센서별 값은 6.4·6.5이 다룬다. 마르코프성은 $p(x(t+\Delta)\mid x(s),s\le t)=p(x(t+\Delta)\mid x(t))$이다. 평균 0, 분산 1로 정규화한 정상 가우시안 과정에서 $\mathbb E[x(t+\tau)\mid x(t)]=\rho(\tau)\,x(t)$이므로 마르코프성은 $\rho(\tau_1+\tau_2)=\rho(\tau_1)\rho(\tau_2)$를 요구하고, 연속인 해는 $\rho(\tau)=e^{-\alpha\lvert\tau\rvert}$이다. Rasmussen·Williams의 부록 B.2.1은 OU 과정이 유일한 정상 1차 가우시안 마르코프 과정이라고 적는다(그 증명은 읽지 않았다 (확인 필요: 출처)). OU 과정의 원 논문은 Uhlenbeck·Ornstein(1930)이고 관련 문헌으로 Doob(1942)이 있다(서지만 확인). OU 과정은 다음 SDE의 정상 해이다.

$$dx=-\alpha x\,dt+\sigma\,dW,\qquad R(\tau)=\frac{\sigma^2}{2\alpha}e^{-\alpha\lvert\tau\rvert},\qquad S(\omega)=\frac{\sigma^2}{\alpha^2+\omega^2}$$

간격 $\Delta t$의 표본은 근사 없이 다음 AR(1)이다. 표본의 PSD는 연속 PSD의 에일리어싱 합 $\sum_kS(\omega+2\pi k/\Delta t)$이며, $\phi=e^{-\alpha\Delta t}$에 대해 $\Delta t\,q/(1-2\phi\cos\omega\Delta t+\phi^2)$와 같다(1.21.5①에서 수치 확인).

$$x_{k+1}=\phi\,x_k+w_k,\qquad \phi=e^{-\alpha\Delta t},\qquad w_k\sim\mathcal N\Bigl(0,\ q\Bigr),\quad q=\frac{\sigma^2}{2\alpha}\bigl(1-\phi^2\bigr)$$

1.21.3.3 가우시안 과정

결론: 가우시안 과정은 평균 함수와 커널로 정해지는 함수의 분포이다. 잡음이 가우시안이면 사후 예측과 주변 우도가 닫힌 형식으로 얻어지며 비용은 학습 점 수의 세제곱이다.

정의는 유한 개의 값이 결합 가우시안인 확률변수의 모음이다(Rasmussen·Williams, 정의 2.1). $f\sim\mathcal{GP}(m,k)$이고 $k(x,x')=\mathbb E[(f(x)-m(x))(f(x')-m(x'))]$는 양반정부호이다. RBF 커널은 $k=\sigma_f^2\exp\bigl(-\tfrac{(x-x')^2}{2\ell^2}\bigr)$이고 Matérn 커널은 $\nu=1/2$에서 $\sigma_f^2e^{-r/\ell}$(OU), $\nu=3/2$에서 $\sigma_f^2\bigl(1+\tfrac{\sqrt3r}{\ell}\bigr)e^{-\sqrt3r/\ell}$이다. 표본 함수가 $k$번 평균제곱 미분 가능할 조건은 $\nu>k$이다(같은 책 4.2). 관측 $y=f(x)+\varepsilon$, $\varepsilon\sim\mathcal N(0,\sigma_n^2)$에서 $K=[k(x_i,x_j)]$, $k_*=[k(x_*,x_i)]$이면 사후 예측은 1.19.3.2의 조건부 공식에 $x_b=y$를 대입한 결과이다.

$$\bar f_*=k_*^\top(K+\sigma_n^2I)^{-1}y,\qquad \mathbb V[f_*]=k(x_*,x_*)-k_*^\top(K+\sigma_n^2I)^{-1}k_*$$

분산은 $y$에 달려 있지 않고 입력 위치로만 정해진다. 관측 $y_*$의 예측 분산은 $\sigma_n^2$을 더한다. 사전 예측은 $f\sim\mathcal N(0,K(X_*,X_*))$이고 사전·사후 표본은 해당 공분산의 촐레스키로 만든다(1.19.3.1). 주변 우도는 $y\sim\mathcal N(0,K+\sigma_n^2I)$에서 얻으며 첫 항이 데이터 적합, 둘째 항이 복잡도 벌점이다(식 2.30, 5.4.1). 하이퍼파라미터 $\theta=(\ell,\sigma_f,\sigma_n)$은 이 값을 최대화해 정한다.

$$\log p(y\mid X,\theta)=-\tfrac12y^\top(K+\sigma_n^2I)^{-1}y-\tfrac12\log\lvert K+\sigma_n^2I\rvert-\tfrac N2\log2\pi$$

계산은 촐레스키 $K+\sigma_n^2I=LL^\top$로 한다(알고리즘 2.1). 비용은 분해에 $N^3/6$, 시험점마다 $N^2/2$이고 기억 용량은 $N^2$이다. RBF 커널 행렬은 점이 가깝거나 $\ell$이 크면 수치적으로 특이하므로 대각에 지터를 더한다. 희소 근사는 $m\ll N$개의 유도점으로 비용을 $O(Nm^2)$로 줄인다. 부분 회귀(SR)는 평균이 $k_m(x_*)^\top(K_{mn}K_{nm}+\sigma_n^2K_{mm})^{-1}K_{mn}y$이고 분산이 $\sigma_n^2k_m^\top(K_{mn}K_{nm}+\sigma_n^2K_{mm})^{-1}k_m$이다(8.3.1, 식 8.14~8.15). 변분 방식(Titsias 2009)과 통합 해석(Quiñonero-Candela·Rasmussen 2005)은 이 장에서 구현하지 않았다. 시간 입력의 Matérn $\nu=p-1/2$ 커널은 연속 시간 AR($p$) 정상 해이다(Rasmussen·Williams 4.2). 이 표현은 칼만 필터·스무더로 $O(N)$에 회귀를 계산하는 방법의 근거이며(Hartikainen·Särkkä 2010; Särkkä·Solin 2019; 서지만 확인) 이 장은 $\nu=1/2$에서 확인했다(1.21.5②). GP 사전으로 연속 시간 궤적을 추정하는 방법은 9.10, 학습 기반 동역학은 12.27, 획득 함수를 사용하는 베이지안 최적화는 1.16가 다룬다.

1.21.3.4 몬테카를로 방법과 중요도 표본추출

결론: 독립 표본 평균의 오차는 차원과 무관하게 $N^{-1/2}$로 줄어든다. 기각 표본추출의 수락률은 차원에 지수적으로 감소하고, 중요도 표본추출은 제안 분포가 목표 분포와 어긋나면 유효 표본 수가 급감한다.

$I=\mathbb E_p[h(x)]$의 표본 평균 $\hat I=\frac1N\sum h(x_i)$는 분산이 $\mathrm{Var}_p[h]/N$이다(1.19.3.5). 이 이름의 방법을 다룬 초기 논문으로 Metropolis·Ulam(1949)이 있다(서지만 확인). 기각 표본추출은 $Mq(x)\ge p(x)$인 제안 $q$에서 $x$를 뽑아 $u\,Mq(x)\le p(x)$일 때 수락하며 수락률이 $1/M$이다. 단위 초구를 정육면체에서 기각 표본추출하면 수락률은 부피비 $\pi^{d/2}/(2^d\Gamma(d/2+1))$이다. 중요도 표본추출은 $\mathbb E_p[h]=\mathbb E_q[hw]$, $w=p/q$를 이용한다. 정규화 상수를 모르는 경우의 자기 정규화 추정량은 $\tilde I=\sum w_ih(x_i)/\sum w_i$이다. 유효 표본 수는 다음과 같다(Owen 식 9.13~9.14; Martino 외 식 6).

$$N_{\mathrm{eff}}=\frac{\bigl(\sum_iw_i\bigr)^2}{\sum_iw_i^2}=\frac{1}{\sum_i\bar w_i^2},\qquad N_{\mathrm{eff}}^*=\frac{N}{\mathbb E_q[w^2]}$$

용어의 사용례로 Kong 외(1994)를 Vehtari 외가 인용한다. 이 식은 추정량 분산 비로 정의한 이론적 유효 표본 수의 근사이며 특정한 경우에만 정확하다(Martino 외). 한 표본 집합의 가중치로 계산하므로 같은 가중치가 만든 문제를 드러내지 못할 수 있다(Owen 9.3). 이 장의 계산에서는 꼬리가 무거울 때 표본 값이 이론 값보다 크게 나왔다. 가우시안 $p=\mathcal N(0,I_d)$, $q=\mathcal N(m\mathbf 1,s^2I_d)$에서는 $s^2>1/2$일 때 $\mathbb E_q[w^2]=\bigl[\tfrac{s^2}{\sqrt{2s^2-1}}e^{m^2/(2s^2-1)}\bigr]^d$이고 그렇지 않으면 무한이다. $m=0$인 경우는 Owen의 예 9.3과 같은 식이다. 따라서 $N_{\mathrm{eff}}/N$은 차원과 평균 이동에 지수적으로 감소한다. 이 현상이 파티클 필터의 가중치 퇴화이며 리샘플링은 9.5, 가중치를 사용하는 표본 기반 제어는 1.16·12.19가 다룬다. 분산 감소는 같은 비용으로 추정량의 분산을 줄이는 기법군이며 대립 변수·제어 변수·층화·준몬테카를로가 있다. 효과는 피적분 함수의 매끄러움과 차원에 달려 있다.

1.21.3.5 마르코프 체인 몬테카를로

결론: 상세 균형을 만족하는 마르코프 체인은 목표 분포를 정상 분포로 가진다. 메트로폴리스-헤이스팅스는 정규화 상수 없이 이를 구현하며 제안 폭이 수락률과 자기상관의 절충을 정한다. R-hat은 사슬 사이의 일치만 검사한다.

현재 값 $x$에서 $x'\sim q(x'\mid x)$를 제안하고 다음 확률로 수락한다(Metropolis 외 1953; Hastings 1970). 대칭 제안이면 $p(x')/p(x)$만 필요하다. 상세 균형 $p(x)T(x\to x')=p(x')T(x'\to x)$가 성립하면 $p$가 정상 분포이다.

$$a(x\to x')=\min\Bigl\{1,\ \frac{p(x')\,q(x\mid x')}{p(x)\,q(x'\mid x)}\Bigr\}$$

랜덤 워크 제안 $x'=x+s\xi$에서 $s$가 작으면 수락률은 높으나 이동이 느리고 $s$가 크면 대부분 기각된다. 독립 성분의 고차원 극한에서 최적 수락률은 약 0.234로 알려져 있다(Roberts 외 1997; (확인 필요: 정리 번호)). 효율은 적분 자기상관 시간 $\tau=1+2\sum_{t\ge1}\rho_t$와 $N_{\mathrm{eff}}=N/\tau$로 잰다(Vehtari 외 식 7). 자기상관의 합은 Geyer의 초기 양의 열로 절단한다(Geyer 1992; arviz-stats 소스의 주석에서 확인). 분할 $\hat R$는 각 사슬을 반으로 나눈 $M$개 조각에서 구간 사이 분산 $B$와 구간 안 분산 $W$로 계산한다(Vehtari 외 식 1~4).

$$\widehat{\mathrm{var}}^+=\frac{N-1}{N}W+\frac1NB,\qquad \hat R=\sqrt{\widehat{\mathrm{var}}^+/W}$$

Vehtari 외는 $\hat R<1.01$과 사슬 4개 이상을 권고하며, 시작 값이 목표 분포보다 넓게 흩어져 있어야 한다고 적는다. $N_{\mathrm{eff}}$는 사후 평균의 추정 효율에만 대응한다는 점도 같은 글이 밝힌다. 따라서 $\hat R$ 값은 효율을 재지 않으며, 모든 사슬이 같은 영역에 머무르면 1에 가깝다.

1.21.4 알고리즘

OU 모의, 가우시안 과정 회귀, 중요도 표본추출, 메트로폴리스-헤이스팅스, 분할 R-hat

ou_exact(α, σ², Δt):  φ ← exp(−α Δt) ;  x₀ ~ N(0, σ²/(2α))
                      x_{k+1} ← φ x_k + √(σ²/(2α)·(1−φ²)) ξ_k
gp_fit(X, y, k):      L ← chol(K + σ_n² I) ;  α ← Lᵀ \ (L \ y)
                      log p(y) ← −½ yᵀα − Σ log L_ii − (N/2) log 2π
gp_pred(x*):          v ← L \ k* ;  평균 ← k*ᵀα ;  분산 ← k(x*,x*) − vᵀv
importance(p, q, h):  x_i ~ q ;  w_i ← p(x_i)/q(x_i) ;  Î ← Σ w_i h(x_i) / Σ w_i
                      N_eff ← (Σ w_i)² / Σ w_i²
metropolis(p, s):     y ← x + s ξ ;  u ~ U(0,1) ;  u < p(y)/p(x)이면 x ← y
split_rhat(chains):   사슬을 반으로 나눠 B, W 계산 ;  √(((N−1)W/N + B/N) / W)

구현은 1.21.5의 코드에 있다. 스펙트럼 추정, 칼만 필터 로그 우도, 희소 근사(SR), 분산 감소의 보조 함수는 스크립트 robotics/fig-src/1.21-process.py에만 있다.

1.21.5 풀이 예제 — OU 과정·가우시안 과정 회귀·중요도 표본추출·메트로폴리스-헤이스팅스와 라이브러리 확인

네 가지를 순수 파이썬으로 계산하고 scikit-learn·SciPy·ArviZ와 같은 입력으로 비교한다. 결론부터 적으면 다음과 같다.

• ① 간격 0.1 s의 OU 표본 $10^5$개 20회에서 표본 자기상관은 이론 $e^{-\alpha\tau}$와 0.0025 이내로 일치한다. 평균의 유효 표본 수는 $N$의 2.5%이며 PSD 추정은 8개 주파수에서 표본화 이론과 1.4% 이내로 일치한다.
• ② 가우시안 과정은 주변 우도 최대에서 RMSE 0.078이고, $\ell=3.0$ 고정에서는 95% 구간이 참 함수의 51%만 포함한다. OU 커널 회귀의 로그 주변 우도는 칼만 필터와 $5.7\times10^{-14}$ 안에서 같다.
• ③ 중요도 표본추출의 이론 $N_{\mathrm{eff}}/N$은 표본 값과 $m\le1$에서 일치하고 $m\ge2$에서는 표본 값이 크게 나온다. 희귀 사건에서는 분산 비 203을 얻는다.
• ④ 상관 가우시안의 $N_{\mathrm{eff}}/N$은 측정한 제안 폭 8개 중 수락률 0.10에서 가장 큰 0.0426이다. 쌍봉 분포에서 네 사슬이 한 봉우리에 머물면 $\hat R=1.000$이다.
• ⑤ scikit-learn과 SciPy는 사후 평균·분산·로그 주변 우도가 $4\times10^{-14}$ 이내로 일치하고, ArviZ의 분할 $\hat R$ 값은 이 장과 같다.

구현. ar1_path·acov·ess_is·is_gauss·mh·iact·split_rhat, rbf·gp_fit·gp_pred. 1.2의 chol·tri_solve·transpose, 1.19의 sample(가우시안 표본), 1.1의 solve_lu(비교)

import math
from operator import mul

def ar1_path(rng, a, s2, dt, n):
    """OU의 정확한 이산화 x' = phi x + w (s2: 정상 분산)"""
    phi = math.exp(-a * dt)
    sq = math.sqrt(s2 * (1 - phi * phi))
    x, out = rng.gauss(0, math.sqrt(s2)), []
    for _ in range(n):
        out.append(x)
        x = phi * x + sq * rng.gauss(0, 1)
    return out

def acov(x, nlag):
    """표본 자기공분산 R[0..nlag] (표본 평균을 뺀다)"""
    m = sum(x) / len(x)
    y = [v - m for v in x]
    return [sum(map(mul, y, y[k:])) / len(y)
            for k in range(nlag + 1)]

def ess_is(w):
    """유효 표본 수 (sum w)^2 / sum w^2 (Owen 식 9.13)"""
    return sum(w) ** 2 / sum(v * v for v in w)

def is_gauss(rng, n, d, m, s):
    """q = N(m, s^2 I_d)의 표본과 p = N(0, I_d)와의 비 p/q"""
    X = [[m + s * rng.gauss(0, 1) for _ in range(d)]
         for _ in range(n)]
    W = [s ** d * math.exp(sum(
        0.5 * ((v - m) / s) ** 2 - 0.5 * v * v for v in x))
         for x in X]
    return X, W

def mh(rng, logp, x0, step, n):
    """랜덤 워크 메트로폴리스(대칭 제안): 수락률과 경로"""
    x, lp, out, acc = list(x0), logp(x0), [], 0
    for _ in range(n):
        y = [v + step * rng.gauss(0, 1) for v in x]
        ly = logp(y)
        if math.log(1.0 - rng.random()) < ly - lp:
            x, lp, acc = y, ly, acc + 1
        out.append(x)
    return out, acc / n

def iact(x, maxlag=3000):
    """적분 자기상관 시간 1 + 2 sum rho_t (Geyer 쌍 합)"""
    m = sum(x) / len(x)
    y = [v - m for v in x]
    r0 = sum(v * v for v in y)
    rho = lambda k: sum(map(mul, y, y[k:])) / r0
    tau, k = -1.0, 0
    while 2 * k + 1 <= maxlag:
        g = rho(2 * k) + rho(2 * k + 1)
        if g <= 0:
            break
        tau, k = tau + 2 * g, k + 1
    return tau

def split_rhat(chains):
    """분할 R-hat: 사슬을 반으로 나눠 B, W로 계산"""
    h = len(chains[0]) // 2
    ch = [c[:h] for c in chains]
    ch += [c[h:2 * h] for c in chains]
    M, N = len(ch), h
    mu = [sum(c) / N for c in ch]
    g = sum(mu) / M
    B = N * sum((v - g) ** 2 for v in mu) / (M - 1)
    W = sum(sum((v - m) ** 2 for v in c) / (N - 1)
            for c, m in zip(ch, mu)) / M
    return math.sqrt(((N - 1) / N * W + B / N) / W)

import math

def rbf(a, b, ell, sf):
    return sf * sf * math.exp(-0.5 * ((a - b) / ell) ** 2)

def gp_fit(X, y, ell, sf, sn):
    """K + sn^2 I = L L^T, alpha, 로그 주변 우도"""
    n = len(X)
    K = [[rbf(a, b, ell, sf) + (sn * sn if i == j else 0.0)
          for j, b in enumerate(X)] for i, a in enumerate(X)]
    L = chol(K)
    al = tri_solve(transpose(L), tri_solve(L, y, True), False)
    lml = (-0.5 * sum(a * b for a, b in zip(y, al))
           - sum(math.log(L[i][i]) for i in range(n))
           - 0.5 * n * math.log(2 * math.pi))
    return L, al, lml

def gp_pred(X, L, al, ell, sf, xs):
    """잠재 함수 f의 사후 평균과 분산 (알고리즘 2.1)"""
    out = []
    for x in xs:
        k = [rbf(x, b, ell, sf) for b in X]
        v = tri_solve(L, k, True)
        out.append((sum(a * b for a, b in zip(k, al)),
                    sf * sf - sum(t * t for t in v)))
    return out

① OU 과정. 상관 시간은 $1/\alpha=2$ s이다. 표본 자기상관은 이론과 다음과 같이 일치한다.

분산 추정값의 평균은 1.0027(이론 1)이고 표준편차는 0.0220(이론 0.0200)이다. 표본 평균의 표준편차 0.0190은 이론 0.0200과 같은 크기이며 화이트 잡음 가정의 $N^{-1/2}=0.0032$보다 6.3배 크다. 평균에 대한 이론 유효 표본 수는 2499로 $N$의 2.5%이다. 위너 과정(20,000 경로)의 앙상블 분산은 $t=1,5,20$ s에서 1.022·4.996·20.330($\sigma^2t$와 비교)으로 선형 증가하고, 같은 잡음 세기의 OU 과정은 0.630·0.990·1.000으로 정상 분산 1에 포화한다. PSD는 자기공분산에 Hann 창 함수를 곱한 Blackman-Tukey 추정(지연 400개)으로 구했다.

추정은 표본화 이론과 일치하고 연속 이론과는 $\omega$가 나이퀴스트 각주파수 $\pi/\Delta t=31.4$ rad/s에 가까울수록 어긋난다. $\omega=30$ rad/s에서 표본화 PSD는 연속 PSD의 2.3배이다. 에일리어싱 합과 닫힌 형식은 400,001항에서 상대 $1.0\times10^{-6}$ 이내로 같다(절단 오차).

② 가우시안 과정 회귀. 하이퍼파라미터를 격자에서 로그 주변 우도 최대로 정하고, 다른 값으로 고정한 설정과 참 함수 대비 RMSE·95% 구간 포함 비율을 비교한다.

$\ell$이 짧으면(0.3) $x\le6$의 평균 사후 표준편차가 0.304로 최대 우도 설정의 0.113보다 크고 RMSE가 0.291이다. $\ell$이 길면(3.0) 데이터를 잡음으로 설명해 $\sigma_n=0.398$이 되고 관측 구간에서 포함 비율이 0.51로 과신한다. $\ell$이 길수록 $\sigma_n$이 커지는 경향은 Rasmussen·Williams 2.3(그림 2.5)에서도 보고된다. 사전 표본 4,000개(격자 21점)의 표본 공분산은 커널 행렬과 최대 0.019 차이가 나고, 사후 표본의 평균·표준편차는 사후 값과 각각 최대 0.0121·0.0232 차이가 난다. 촐레스키는 RBF 커널 행렬($N=50$, $\ell=1$)에서 지터 $10^{-15}$까지 성공하고 $10^{-16}$에서 실패한다. 시간 입력의 OU 커널은 칼만 필터와 같은 로그 주변 우도를 주고 비용이 선형이다.

촐레스키의 CPU 시간은 $N$이 2배가 될 때 7.1·8.0·9.3배로 늘어($N=100\to800$: 0.011, 0.080, 0.639, 5.926 s) 이론 8에 가깝다. 시간은 실행마다 달라지며 이 값은 출력 파일의 실행 하나이다. 같은 $N=400$ 자료에 SR을 적용하면 다음과 같다.

유도점 40개의 SR은 사후 평균이 전체 GP와 0.0006 이내이고 CPU 시간이 약 1/40(0.04 s 대 1.50 s)이다. 유도점에서 먼 $x=8$에서는 사후 표준편차가 0.538로 전체 GP(0.955)의 56%이며 $m=5$에서는 0.001이다. SR의 분산 공식이 유도점에서 멀수록 0으로 가는 구조 때문이다.

③ 몬테카를로. 정육면체에서 초구를 기각 표본추출하면 수락률은 차원 2·5·10에서 모의 0.786·0.162·0.0026(이론 0.785·0.165·0.0025)이고 $d=20$에서는 이론 $2.5\times10^{-8}$이며 $10^5$회 시도에서 수락 0회이다. 중요도 표본추출에서 $\mathbb E[x^2]=1$을 자기 정규화 추정한다.

표본 $N_{\mathrm{eff}}/N$은 가중치가 균일에 가까울 때($m\le1$과 $s=1.5$) 이론과 일치하고 꼬리가 무거워지면 이론보다 크다. $m=3$의 표본 $N_{\mathrm{eff}}=42$로 예측한 RMSE는 $\sqrt{2/42}=0.22$이나 실제는 0.478이다. $s=0.5$는 $\mathbb E_q[w^2]$가 무한한데도 표본 $N_{\mathrm{eff}}/N$이 0.0992로 나온다. 차원에 따른 퇴화는 다음과 같다.

희귀 사건 $P(X>3)=1.3499\times10^{-3}$을 $N=10^4$로 200회 추정하면 단순 몬테카를로의 상대 표준편차는 0.268, $q=\mathcal N(3,1)$의 중요도 표본추출은 0.019로 분산 비가 203이다. 이 제안의 전체 표본 $N_{\mathrm{eff}}/N$은 0.0041에 불과하다. 따라서 $N_{\mathrm{eff}}$는 피적분 함수와 무관한 지표이고 이 추정량의 효율을 대표하지 못한다. Owen도 중요도 표본추출의 효율이 $f$에 달려 있다고 적는다. 분산 감소 기법은 $\int_0^1e^x\,dx$에서 비교했다(단순 몬테카를로의 이론 표준편차 $1.556\times10^{-2}$).

④ 메트로폴리스-헤이스팅스. $\rho=0.9$인 2차원 가우시안에서 제안 폭을 바꾼다. $N_{\mathrm{eff}}$는 첫 사슬의 $x_1$에서 구했다(워밍업 2,000개 제외).

$N_{\mathrm{eff}}/N$은 수락률 0.10~0.31($s=1.0\sim2.5$)에서 0.033~0.043으로 가장 크고, 수락률이 가장 높은 $s=0.1$에서는 0.0012(사슬당 $N_{\mathrm{eff}}$ 약 46)이다. 이때도 분할 $\hat R$는 1.006으로 임계 1.01 아래이므로 $\hat R$ 값은 효율을 재지 않는다. 쌍봉 목표(간격 12)에서 수렴 진단의 한계를 본다.

가에서는 네 사슬이 모두 오른쪽 봉우리에 머물러 $\hat R=1.000$이고 사슬 하나의 $N_{\mathrm{eff}}$ 1768도 양호하나 평균이 +5.99로 봉우리 하나가 누락된다. 나에서는 $\hat R=6.48$이 이를 검출하지만 두 사슬씩이라 평균이 +0.02, 비율이 0.50으로 우연히 참값에 가깝다. 다에서는 폭 8의 제안이 봉우리 사이를 사슬당 194.8회 이동하고 $\hat R=1.007$이다. 사슬 하나의 $N_{\mathrm{eff}}$는 가·나에서 1768·1821로 양호해 두 경우를 구분하지 못한다.

⑤ 라이브러리 호출 결과(실행 2026-11-03). 스크립트는 robotics/fig-src/1.21-libcheck.py, 출력은 robotics/fig-src/data/1.21-libcheck-output.txt에 있다.

WhiteKernel을 사용하면 return_std가 관측 잡음을 포함한 예측 표준편차를 돌려준다. alpha를 사용하는 방식에서는 잠재 함수의 표준편차이다. 소스에서 predict가 kernel_.diag(X)를 사용하고 WhiteKernel.diag가 잡음 수준을 돌려주는 것을 읽었다. 두 방식의 평균과 우도는 같다. scikit-learn의 최적화는 재시작 없이도 재시작 10회와 같은 해에 도달했다. 순위 정규화 $\hat R$는 가에서 1.0003으로 이 경우도 검출하지 못하고, 나에서는 6.48보다 작은 1.73이나 1.01을 넘는다. 다중 사슬 ESS가 나에서 6으로 봉우리 수에 가깝다는 점은 Vehtari 외의 서술과 부합한다. numpy의 단일 스레드 촐레스키는 $N=1000,2000,4000$에서 0.026, 0.129, 1.035 s이다.

한계. 수치는 합성 입력(시드 고정)이다. GP는 한 데이터 집합과 RBF 커널만 다뤘고 하이퍼파라미터는 격자(로그 간격 약 20%)로 탐색했으며 포함 비율은 참 함수를 아는 경우에 계산했다. 중요도 표본추출은 가우시안 쌍(해석값이 있는 경우)만, MH는 랜덤 워크 제안과 2차원 사례만 다뤘다. 쌍봉 사례는 진단의 한계를 보이도록 구성한 것이다. 시간은 부하 평균 17~20의 공유 기계에서 측정한 CPU 시간(최소값)이다.

견해: 이 책은 가우시안 과정의 하이퍼파라미터를 주변 우도 최대화로 정한 뒤 보류 관측의 예측 구간 포함 비율로 점검하는 것을 기본으로 둔다고 본다. 근거는 $\ell=3.0$ 고정 설정이 로그 주변 우도가 최대보다 6.98 낮고 95% 구간 포함 비율이 0.51로 과신이 드러났다는 점이다. 참 함수를 알 수 없는 실제 데이터에서는 보류 관측으로 대신해야 하며 이 장은 참 함수로 점검했다. $N=12$의 한 사례이므로 우도 곡면이 평평한 경우는 측정하지 않았다 (추정).

1.21.6 설계 트레이드오프 — 불확실성 모델과 계산 방식

공개 문헌에서 확인되는 선택지는 네 가지다.

• ① 시간 상관 잡음을 화이트 노이즈로 근사. 구현이 가장 단순하다. 상관 시간이 표본 간격보다 길면 평균의 정확도를 과대 평가한다.
• ② 1차 가우스-마르코프(OU) 상태를 증강. 상태 하나를 추가하며 칼만 필터에서 선형 비용이다(9.3). 상관 시간과 분산 두 모수를 식별해야 한다.
• ③ 가우시안 과정 모델. 함수 형태를 가정하지 않고 불확실성을 준다. 비용은 $O(N^3)$이며 희소·상태공간 근사가 필요하다(9.10, 12.27).
• ④ 표본 기반(중요도 표본추출·MCMC). 임의 분포를 다룬다. 가중치 퇴화와 수렴 진단이 필요하다(9.5).

견해: 시간 상관이 확인된 센서 잡음은 화이트 노이즈로 두지 않고 유효 표본 수나 OU 상태로 반영하는 것이 맞다고 본다. 근거는 이 장의 OU 예에서 $N=10^5$개 평균의 표준편차가 0.0200으로 화이트 가정의 0.0032보다 6.3배 크다는 점이다. 상관 시간 2 s·표본 간격 0.1 s라는 한 조건이며 다른 비율에서는 배율이 달라진다.

1.21.7 응용 축 비교

이 표는 이 책의 정리이며 장별 요구는 다르다.

1.21.8 관련 규격·프로토콜

통계 용어와 기호에 대한 ISO 3534 계열이 있는 것으로 알려져 있으나 이 장은 본문을 확인하지 않았다 (확인 필요). 알고리즘의 상호운용 기준은 정의와 하이퍼파라미터 규약(예: scikit-learn의 alpha와 WhiteKernel의 차이)을 사용하는 쪽의 문서에서 읽는 것이다.

1.21.9 관련 규제·인증 요구와의 연결

직접 적용되는 규제는 없다. 몬테카를로 시험·시뮬레이션 결과를 안전 논증의 근거로 제출할 때는 표본 수와 신뢰 구간이 심사 항목이 될 수 있으며 통계적 안전 입증은 26.11, 시나리오 기반 시험은 16.24, 시뮬레이션 신뢰성 입증은 41.14이 다룬다. 마르코프 체인의 수렴을 R-hat만으로 판단하는 것은 확률적 보증을 내세우는 소프트웨어를 검증할 때 점검 항목이 될 수 있다.

1.21.10 한계·유의점·미해결

• 정리의 출처. 위너-힌친 정리, 위너 과정의 비미분성, OU 과정의 유일성, 최적 수락률 0.234의 정리 번호와 교과서 장 번호는 확인하지 않았다 (확인 필요).
• 가우시안 과정. RBF 한 종류와 OU(Matérn $\nu=1/2$)만 계산했다. 다차원 입력·ARD·비정상 커널·분류는 다루지 않았다. 희소 근사는 SR만 구현했고 변분 방식은 구현하지 않았다.
• 몬테카를로. 분산 감소는 1차원 적분 하나로 비교했다. 가중치 꼬리 진단(파레토 평활화 등)은 계산하지 않았다. 순차 몬테카를로는 9.5에 둔다.
• MCMC. 랜덤 워크 제안만 다뤘고 해밀턴 몬테카를로·깁스·적응형 제안은 다루지 않았다.
• 라이브러리. GPy·GPyTorch·emcee·Stan은 실행하지 않았다. scikit-learn과 arviz-stats는 소스의 일부만 읽었다.

1.21.11 체크포인트·연구 질문

체크포인트:

1. 광의 정상 과정의 자기공분산에서 PSD를 구하고, OU 과정의 두 식을 유도할 수 있는가.
2. 상관된 표본 평균의 유효 표본 수를 자기상관으로 계산하고 화이트 노이즈 가정의 오차 배율을 말할 수 있는가.
3. 가우시안 과정 사후 평균·분산·주변 우도를 촐레스키로 적을 수 있는가.
4. $\ell$이 너무 짧거나 길 때 구간이 어떻게 달라지는지, 지터가 필요한 이유를 설명할 수 있는가.
5. 중요도 표본추출의 $N_{\mathrm{eff}}$ 식과 그 한계(꼬리, 함수 의존성)를 설명할 수 있는가.
6. 메트로폴리스-헤이스팅스의 수락 확률을 유도하고 $\hat R$가 놓치는 경우를 말할 수 있는가.

연구 질문(47.12에 수록한다):

• 로봇 센서 잡음에서 1차 가우스-마르코프 모델로 충분한 구간은 어디이며 Allan 분산이나 자기상관의 어떤 지표로 판정하는가.
• 가우시안 과정 동역학 모델의 하이퍼파라미터 불확실성은 확률적 MPC의 기회 제약(12.18) 만족에 얼마나 영향을 주는가.
• 중요도 표본추출 기반 제어(12.19)의 유효 표본 수를 운용 중에 어떤 진단으로 감시할 수 있는가.

1차 출처

• C. E. Rasmussen, C. K. I. Williams, 『Gaussian Processes for Machine Learning』, MIT Press, 2006(Crossref 등록 2005, doi:10.7551/mitpress/3206.001.0001) — 공개 PDF의 위 범위를 읽었다.
• A. Vehtari, A. Gelman, D. Simpson, B. Carpenter, P.-C. Bürkner, 「Rank-Normalization, Folding, and Localization: An Improved $\hat R$ for Assessing Convergence of MCMC」, Bayesian Analysis 16(2), 667–718, 2021(doi:10.1214/20-ba1221; arXiv:1903.08008v5) — §1, §3.1~3.2 일부, 권고 부분을 읽었다.
• L. Martino, V. Elvira, F. Louzada, 「Effective sample size for importance sampling based on discrepancy measures」, Signal Processing 131, 386–401, 2017(doi:10.1016/j.sigpro.2016.08.025; arXiv:1602.03572v5) — 초록과 §1~3을 읽었다.
• A. Kong, J. S. Liu, W. H. Wong, 「Sequential Imputations and Bayesian Missing Data Problems」, JASA 89(425), 278–288, 1994(doi:10.1080/01621459.1994.10476469) — 서지만 확인(Vehtari 외가 유효 표본 수 용어의 사용례로 인용).
• A. B. Owen, 『Monte Carlo theory, methods and examples』(온라인 초고, https://artowen.su.domains/mc/) 9장 — 9.2~9.3과 예 9.3을 읽었다.
• JCGM 101:2008, 「Evaluation of measurement data — Supplement 1 to the "Guide to the expression of uncertainty in measurement" — Propagation of distributions using a Monte Carlo method」 — §1, §7.2를 읽었다.
• N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth, A. H. Teller, E. Teller, 「Equation of State Calculations by Fast Computing Machines」, J. Chem. Phys. 21(6), 1087–1092, 1953(doi:10.1063/1.1699114); W. K. Hastings, Biometrika 57(1), 97–109, 1970(doi:10.1093/biomet/57.1.97); N. Metropolis, S. Ulam, JASA 44(247), 335–341, 1949(doi:10.1080/01621459.1949.10483310) — 서지만 확인.
• A. Gelman, D. B. Rubin, Statistical Science 7(4), 457–472, 1992(doi:10.1214/ss/1177011136); C. J. Geyer, Statistical Science 7(4), 473–483, 1992(doi:10.1214/ss/1177011137); G. O. Roberts, A. Gelman, W. R. Gilks, Annals of Applied Probability 7(1), 110–120, 1997(doi:10.1214/aoap/1034625254) — 서지만 확인.
• G. E. Uhlenbeck, L. S. Ornstein, 「On the Theory of the Brownian Motion」, Phys. Rev. 36(5), 823–841, 1930(doi:10.1103/physrev.36.823); J. L. Doob, 「The Brownian Movement and Stochastic Equations」, Annals of Mathematics 43(2), 1942(doi:10.2307/1968873); N. Wiener, 「Generalized harmonic analysis」, Acta Mathematica 55, 117–258, 1930(doi:10.1007/bf02546511); A. Khintchine, Mathematische Annalen 109(1), 604–615, 1934(doi:10.1007/bf01449156) — 서지만 확인.
• J. Quiñonero-Candela, C. E. Rasmussen, 「A Unifying View of Sparse Approximate Gaussian Process Regression」, JMLR 6, 1939–1959, 2005; M. Titsias, 「Variational Learning of Inducing Variables in Sparse Gaussian Processes」, PMLR 5(AISTATS), 567–574, 2009; J. Hartikainen, S. Särkkä, MLSP 2010, 379–384(doi:10.1109/mlsp.2010.5589113); S. Särkkä, A. Solin, 『Applied Stochastic Differential Equations』, Cambridge University Press, 2019(doi:10.1017/9781108186735) — 서지만 확인(JMLR·PMLR은 출판사 페이지).
• scikit-learn 1.9.1(GaussianProcessRegressor·Matern·WhiteKernel), SciPy 1.18.1(stats·signal.welch), NumPy 2.5.4, ArviZ 1.3.0·arviz-stats 1.3.3(rhat·ess) 실행, 2026-11-03.

다음 장. 1.22은 구간 해석과 집합 기반 추정으로, 구간 연산·수축자·보장된 위치추정을 다룬다. 이 장까지의 확률 모델은 오차를 분포로 둔다. 1.22은 오차가 경계 안에 있다는 가정만으로 위치를 구간 또는 집합으로 보장하는 방식을 정리한다.