E10 · 모의 경로를 뽑는 일과 베이지안 추론의 선형대수

내용

E10 · 모의 경로를 뽑는 일과 베이지안 추론의 선형대수#

두 산업의 충격을 함께 뽑으려면#

두 산업의 다음 분기 생산 증가율에서 이미 예측한 평균을 뺀 값을 $x_1,x_2$라 합시다. 단위는 퍼센트포인트이고, 공분산의 단위는 그 제곱입니다. 이번 장에서는 먼저

\[\begin{split} x\sim N(0,\Sigma),\qquad \Sigma=\begin{pmatrix}4&2\\2&3\end{pmatrix} \end{split}\]

를 가정합니다. 두 개의 독립 표준정규 난수를 그대로 쓰면 공분산은 $I$입니다. 원하는 충격을 만들려면 독립 난수에 서로 다른 크기와 공통 성분을 주어야 합니다. $L=\left(\begin{smallmatrix}a&0\\b&c\end{smallmatrix}\right)$라 놓고 $LL^T=\Sigma$의 세 성분을 맞추면

\[ a^2=4,\quad ab=2,\quad b^2+c^2=3. \]

양의 대각을 선택하면 $a=2,b=1,c=\sqrt2$입니다. 따라서

\[ x_1=2z_1,\qquad x_2=z_1+\sqrt2z_2,\qquad z_1,z_2\overset{\rm iid}{\sim}N(0,1). \]

첫 충격의 절반을 둘째 산업이 공유합니다. $z=(1,-1)^T$가 뽑힌 한 번의 실험에서는 $x=(2,1-\sqrt2)^T$입니다. 한 번의 표본이 공분산을 재현해야 하는 것은 아닙니다. 반복해서 뽑은 분포의 공분산이 $LL^T$라는 뜻입니다.

같은 축에 그린 원형 표준정규 등고선과 공분산 4 2 2 3의 타원, 그리고 특이 공분산의 직선 지지집합 — 그림 158 반지름 1인 원을 $L$로 보낸 타원이다. 확률질량 전체의 경계가 아니라 Mahalanobis 거리 1인 등고선이다. 점선은 뒤에서 다룰 rank 1 공분산의 지지 직선이다. 두 축의 단위와 눈금은 같다.#

이 표본생성법의 조건은 $AA^T=\Sigma$뿐입니다. 고유분해 $U\Lambda U^T$에서 $A=U\Lambda^{1/2}$를 써도 됩니다. 양정치 조밀행렬에서는 Cholesky가 대개 더 적은 연산으로 끝나며, 분해에 $O(d^3)$, 표본마다 삼각행렬 곱에 $O(d^2)$가 듭니다. 고유분해가 틀린 방법은 아닙니다. rank 판정과 작은 고윳값의 진단에는 오히려 유용합니다. 구현·행렬 구조를 지정하지 않고 고정된 실행시간 배율을 주장하지 않습니다.

$\Sigma_s=\left(\begin{smallmatrix}4&2\\2&1\end{smallmatrix}\right)$이면 $x=(2z,z)$입니다. 평면 위 밀도는 없고 $x_1=2x_2$ 위에 분포합니다. 완전 피벗 Cholesky는 잔여 대각이 가장 큰 좌표를 먼저 선택하고, 양의 피벗만큼 인수를 만듭니다. 정확한 PSD 입력에서 잔여 대각이 모두 0이면 잔여행렬도 0입니다. 실제로 PSD 행렬은 $|r_{ij}|^2\le r_{ii}r_{jj}$이므로 그렇습니다. 수치적으로는 허용오차와 남은 잔차를 함께 보고합니다. 작은 음의 고윳값을 0으로 바꾸는 것은 수정된 공분산에서 추출하는 일이며, 원래 입력의 정확한 표본추출이라고 부르지 않습니다.

관측을 조건으로 넣으면 어느 행렬이 간단해지는가#

두 번째 산업의 증가율이 $x_2=3$으로 관측되었다고 하자. 공분산에서 읽는 조건부 평균과 분산은

\[ E[x_1\mid x_2=3]=\frac23\,3=2,\qquad \operatorname{Var}(x_1\mid x_2)=4-\frac{2^2}{3}=\frac83. \]

이제 정밀도 $\Omega=\Sigma^{-1}=\frac18\left(\begin{smallmatrix}3&-2\\-2&4\end{smallmatrix}\right)$에서 같은 답을 읽어 보겠습니다. $x_2$를 고정한 밀도의 지수는

\[ -\frac12\left(\frac38x_1^2-\frac12x_1x_2\right) =-\frac12\frac38\left(x_1-\frac23x_2\right)^2+\text{$x_1$과 무관한 항}. \]

따라서 조건부 분산은 $\Omega_{11}^{-1}=8/3$입니다. 주변분산은 $\Sigma$의 대각 블록을 읽고, 조건부분산은 $\Omega$의 대각 블록을 역으로 읽습니다. 주변화된 정밀도 자체는 $\Omega_{11}-\Omega_{12}\Omega_{22}^{-1}\Omega_{21}$입니다. 주변화와 조건화를 서로 바꾸지 말자.

평균이 $\mu$인 양정치 Gaussian을 $a,b$ 블록으로 나누면

\[ \mu_{a\mid b}=\mu_a-\Omega_{aa}^{-1}\Omega_{ab}(x_b-\mu_b), \qquad \Sigma_{a\mid b}=\Omega_{aa}^{-1} =\Sigma_{aa}-\Sigma_{ab}\Sigma_{bb}^{-1}\Sigma_{ba}. \]

이 항등식은 H7의 Schur 보원의 통계적 해석입니다. 특이 Gaussian에서는 역행렬 대신 지지공간과 조건부 분포의 존재 범위를 먼저 정해야 합니다. 위 양정치 공식에 임의로 0의 역수를 넣지 않습니다.

세 변수의 예로

\[\begin{split} \Omega=\begin{pmatrix}2&-1&0\\-1&2&-1\\0&-1&2\end{pmatrix}, \qquad \Sigma=\frac14\begin{pmatrix}3&2&1\\2&4&2\\1&2&3\end{pmatrix} \end{split}\]

를 보겠습니다. $x_1,x_3$의 주변공분산은 $1/4$이지만 $x_2$를 알면 조건부 공분산은 $1/4-(1/2)(1)(1/2)=0$입니다. 양정치 Gaussian에서는 $\Omega_{ij}=0$이 나머지 모든 좌표를 조건으로 한 독립과 동치입니다. 일반 분포에서 무상관이 독립을 뜻한다는 주장은 하지 않습니다.

시계열 전체를 한 번에 뽑기#

숨은 경기 수준 $x_t$와 발표 통계 $y_t$를 같은 표준화 단위로 측정하자. 가장 작은 예는

\[ x_1\sim N(0,1),\quad x_{t+1}=x_t+w_t,\quad y_t=x_t+v_t, \qquad w_t,v_t\overset{\rm independent}{\sim}N(0,1). \]

초기 상태와 모든 잡음도 독립입니다. $T=3,y=(1,0,1)^T$일 때 사후 음의 로그밀도의 두 배에서 상수를 버리면

\[ x_1^2+(x_2-x_1)^2+(x_3-x_2)^2+ (1-x_1)^2+x_2^2+(1-x_3)^2 =x^T\Omega x-2h^Tx+2, \]

\[\begin{split} \Omega=\begin{pmatrix}3&-1&0\\-1&3&-1\\0&-1&2\end{pmatrix},\quad h=(1,0,1)^T,\quad m=\Omega^{-1}h=\frac1{13}(6,5,9)^T. \end{split}\]

검산하면 $3(6)-5=13$, $-6+3(5)-9=0$, $-5+2(9)=13$입니다. 사후 공분산은

\[\begin{split} C=\frac1{13}\begin{pmatrix}5&2&1\\2&6&3\\1&3&8\end{pmatrix}. \end{split}\]

서로 이웃한 상태만 동학 방정식에 함께 들어가므로 정밀도는 삼중대각입니다. 공분산은 조밀해도 정밀도의 0은 그대로 남습니다.

세 시점의 관측 1 0 1과 사후 평균 6/13 5/13 9/13, 사후 표준편차 구간을 비교한 그림 — 그림 159 사각형은 발표 통계, 원은 세 관측을 모두 사용한 사후 평균이다. 세로선은 각 좌표의 사후 평균 ±1 표준편차이며 동시 신뢰띠가 아니다. 이산 시점의 선은 순서를 보여 주기 위한 연결이다.#

스칼라 삼중대각 정밀도의 대각을 $a_t$, 아래대각을 $b_t$라 하면 $\Omega=LDL^T$를

\[ d_1=a_1,\quad \ell_t=b_{t-1}/d_{t-1},\quad d_t=a_t-\ell_t^2d_{t-1} \]

로 구합니다. 예에서는 $d=(3,8/3,13/8)$, $\ell=(-1/3,-3/8)$이고 $\det\Omega=13$입니다. $m$은 앞·뒤 대입으로 구합니다. 독립 $z\sim N(0,I)$에 대해

\[ L^T\eta=D^{-1/2}z,\qquad x=m+\eta \]

로 뽑으면 $\operatorname{Cov}(\eta)=L^{-T}D^{-1}L^{-1}=\Omega^{-1}$입니다. $Lz$를 곱하는 공분산 추출과 달리, 정밀도에서는 전치 삼각계를 푼다.

일반 상태차원 $d$에서 $x_1\sim N(m_0,P_0)$, $x_{t+1}=A_tx_t+w_t$, $y_t=H_tx_t+v_t$이며 $P_0,Q_t,R_t\succ0$라 합시다. 처음·내부·마지막 대각 블록은 각각

\[\begin{split} \begin{aligned} D_1&=P_0^{-1}+A_1^TQ_1^{-1}A_1+H_1^TR_1^{-1}H_1,\\ D_t&=Q_{t-1}^{-1}+A_t^TQ_t^{-1}A_t+H_t^TR_t^{-1}H_t,\\ D_T&=Q_{T-1}^{-1}+H_T^TR_T^{-1}H_T, \end{aligned} \end{split}\]

아래대각은 $-Q_t^{-1}A_t$입니다. $T=1$에는 동학항이 없으므로 $D_1=P_0^{-1}+H_1^TR_1^{-1}H_1$만 씁니다. 정보벡터는 $h_t=H_t^TR_t^{-1}y_t$이며 첫 항에 $P_0^{-1}m_0$를 더합니다. 블록 소거의 비용은 $O(Td^3)$, 저장량은 $O(Td^2)$입니다. $O(T)$라는 말은 $d$가 고정되었다는 뜻입니다.

E8의 Gaussian 필터가 준 $m_t,P_t$와 다음 시점 예측 $a_{t+1},R_{t+1}$로도 뽑을 수 있습니다. 마지막 $x_T\mid y_{1:T}$부터 뽑고

\[ J_t=P_tA_t^TR_{t+1}^{-1},\qquad x_t\mid x_{t+1},y_{1:T}\sim N\bigl(m_t+J_t(x_{t+1}-a_{t+1}),P_t-J_tR_{t+1}J_t^T\bigr) \]

를 뒤로 적용합니다. 조건부 밀도를 연쇄곱하면 같은 사후 결합밀도가 됩니다. 따라서 두 방법은 표본 하나하나가 아니라 분포가 같습니다. 아래 코드는 큰 Monte Carlo 오차에 기대지 않고 두 방법의 평균과 공분산을 직접 비교합니다.

import numpy as np
from scipy.linalg import solve_triangular
O = np.array([[3.,-1,0],[-1,3,-1],[0,-1,2]])
h = np.array([1.,0,1]); y = h.copy()
L = np.linalg.cholesky(O)
m = np.linalg.solve(O,h)
B = solve_triangular(L.T,np.eye(3),lower=False)
C = B@B.T
mf=[]; pf=[]; pred=[]; a=0.; r=1.
for obs in y:
    pred.append(r); k=r/(r+1)
    a=a+k*(obs-a); r=(1-k)*r
    mf.append(a); pf.append(r); r=r+1
ms=np.array(mf); cs=np.zeros((3,3)); cs[-1,-1]=pf[-1]
for t in [1,0]:
    j=pf[t]/pred[t+1]
    ms[t]=mf[t]+j*(ms[t+1]-mf[t])
    cs[t,t]=pf[t]+j*j*(cs[t+1,t+1]-pred[t+1])
    cs[t,t+1:]=j*cs[t+1,t+1:]; cs[t+1:,t]=cs[t,t+1:]
assert np.allclose(m,[6/13,5/13,9/13])
assert np.allclose(ms,m) and np.allclose(cs,C)
rng=np.random.default_rng(410)
draws=m[:,None]+B@rng.standard_normal((3,50000))
print('mean:',m,'; max analytic covariance discrepancy:',abs(cs-C).max())
print('Monte Carlo mean error:',abs(draws.mean(axis=1)-m).max())

# Scalar tridiagonal precision: O(T) storage and arithmetic.
def factor_path(diagonal, offdiagonal):
    d=np.array(diagonal,dtype=float,copy=True)
    ell=np.empty(len(d)-1)
    for t in range(1,len(d)):
        if d[t-1]<=0: raise ValueError('precision is not positive definite')
        ell[t-1]=offdiagonal[t-1]/d[t-1]
        d[t]-=ell[t-1]*offdiagonal[t-1]
    if d[-1]<=0: raise ValueError('precision is not positive definite')
    return d,ell
def backward_path(rhs,ell):
    out=np.array(rhs,dtype=float,copy=True)
    for t in range(len(out)-2,-1,-1): out[t]-=ell[t]*out[t+1]
    return out
def mean_path(h,d,ell):
    z=np.array(h,dtype=float,copy=True)
    for t in range(1,len(z)): z[t]-=ell[t-1]*z[t-1]
    return backward_path(z/d,ell)
d,ell=factor_path([3,3,2],[-1,-1])
assert np.allclose(d,[3,8/3,13/8])
assert np.allclose(mean_path(h,d,ell),m)
# Form the small sampling factor to verify its covariance exactly.
Bpath=np.column_stack([backward_path(np.eye(3)[:,j]/np.sqrt(d),ell) for j in range(3)])
assert np.allclose(Bpath@Bpath.T,C)
T=10000; diagonal=np.r_[np.full(T-1,3.),2.]; off=-np.ones(T-1)
d,ell=factor_path(diagonal,off); observations=rng.normal(size=T)
longmean=mean_path(observations,d,ell)
res=diagonal*longmean-observations
res[1:]+=off*longmean[:-1]; res[:-1]+=off*longmean[1:]
assert np.max(abs(res))<1e-12
pathdraw=longmean+backward_path(rng.normal(size=T)/np.sqrt(d),ell)
assert np.all(np.isfinite(pathdraw))
print('T=10000 precision residual:',np.max(abs(res)))

mean: [0.46153846 0.38461538 0.69230769] ; max analytic covariance discrepancy: 1.1102230246251565e-16
Monte Carlo mean error: 0.005133649128292894
T=10000 precision residual: 1.1102230246251565e-15

베이지안 회귀는 정밀도에 정보를 더한다#

$\beta\sim N(b_0,V_0)$이고 $y\mid\beta\sim N(X\beta,R)$, $V_0,R\succ0$라 합시다. 사전과 우도의 제곱을 전개하면

\[ K=V_0^{-1}+X^TR^{-1}X,\quad h=V_0^{-1}b_0+X^TR^{-1}y, \quad V=K^{-1},\quad b=Vh. \]

한 관측 $X=(1,1)$, $V_0=I$, $R=1$, $b_0=0$, $y=3$이면

\[\begin{split} K=\begin{pmatrix}2&1\\1&2\end{pmatrix},\quad V=\frac13\begin{pmatrix}2&-1\\-1&2\end{pmatrix},\quad b=(1,1)^T. \end{split}\]

관측이 계수의 합을 알려 주기 때문에 두 계수의 사후 오차는 음의 상관을 갖습니다. 관측이 알려 주지 않은 차이 방향의 사전분산은 남습니다. 같은 결과를 관측공간에서 쓰면

\[\begin{split} \begin{aligned} V&=V_0-V_0X^T(R+XV_0X^T)^{-1}XV_0,\\ b&=b_0+V_0X^T(R+XV_0X^T)^{-1}(y-Xb_0). \end{aligned} \end{split}\]

$k$개 계수와 $n\ll k$개 관측일 때 $n\times n$ 계를 풀 수 있다는 이점이 있습니다. 다만 $O(n^3+n^2k)$라는 비용은 $V_0$가 대각이거나 빠르게 곱해지는 경우의 평균·관련 계산에 해당합니다. 조밀한 $V_0$와 전체 $k\times k$ 사후 공분산을 명시적으로 출력하면 그 곱셈·저장 비용도 추가됩니다.

여러 방정식, Wishart, 질량행렬#

$Y=XB+E$에서 $Y$는 $n\times m$, $X$는 $n\times k$입니다. 열 단위 vec 규약으로 $\operatorname{Cov}(\operatorname{vec}E)=\Sigma\otimes I_n$라 합시다. 켤레 사전은

\[ B\mid\Sigma\sim MN(B_0,V_0,\Sigma),\qquad \Sigma\sim IW(S_0,\nu_0),\quad \nu_0>m-1. \]

여기서 역 Wishart 밀도는 $|\Sigma|^{-(\nu_0+m+1)/2}\exp[-\operatorname{tr}(S_0\Sigma^{-1})/2]$에 비례하고 $S_0\succ0$입니다. 행렬 완전제곱은

\[\begin{split} \begin{aligned} V_n^{-1}&=V_0^{-1}+X^TX,& B_n&=V_n(V_0^{-1}B_0+X^TY),\\ S_n&=S_0+(Y-XB_n)^T(Y-XB_n)\\ &\quad +(B_n-B_0)^TV_0^{-1}(B_n-B_0),&\nu_n&=\nu_0+n. \end{aligned} \end{split}\]

를 줍니다. $B$의 조건부 공분산은 $\Sigma\otimes V_n$입니다. $L_VL_V^T=V_n,L_\Sigma L_\Sigma^T=\Sigma$이면 독립 표준정규 행렬 $Z$로 $B=B_n+L_VZL_\Sigma^T$를 뽑습니다. 정규화 상수를 적분할 때 $B-B_n=L_VZL_\Sigma^T$의 Jacobian $|V_n|^{m/2}|\Sigma|^{k/2}$가 $|\Sigma|^{-k/2}$를 상쇄하므로 자유도 갱신은 $\nu_0+n$입니다.

Wishart 계산에 등장하는 두 Jacobian도 확인해 봅시다. 양의 대각 하삼각 $L$에서 $W=LL^T$로 가는 변환은

\[ dW=2^m\prod_{i=1}^m l_{ii}^{m+1-i}\,dL. \]

마지막 행 $v^T,l$를 붙이면 새 비대각은 $L_{m-1}v$, 마지막 대각은 $v^Tv+l^2$입니다. 이전 성분을 먼저 놓은 Jacobian의 새 블록 행렬식은 $\det L_{m-1}\,2l$이고, 귀납적으로 위 곱이 나옵니다.

대칭행렬 위 합동 $H\mapsto AHA^T$의 Jacobian은 $|\det A|^{m+1}$입니다. 대각 $A$에서 대각 좌표의 인자 $a_i^2$와 비대각의 $a_ia_j$를 곱하면 되며, 일반 $A$는 SVD와 직교변환의 부피 보존으로 따릅니다. 역변환 $\Sigma\mapsto\Sigma^{-1}$의 미분 $H\mapsto-\Sigma^{-1}H\Sigma^{-1}$은 따라서 절댓값 Jacobian $|\Sigma|^{-(m+1)}$을 갖습니다. Wishart의 Gaussian 표본 기원과 정규화 상수는 다변량 확률론의 출발 전제로 두고, 여기서는 갱신에 쓰이는 변환을 증명했습니다.

BVAR의 공통 행 공분산 구조를 유지하면 인수분해는 $O(k^3+m^3)$, 표본 곱은 $O(k^2m+km^2)$입니다. 임의의 개별 사전분산은 이 Kronecker 구조를 깨뜨릴 수 있습니다. 공통 더미 관측 $X_0,Y_0$는 $X_0^TX_0=V_0^{-1}$와 $X_0^TY_0=V_0^{-1}B_0$를 구현하지만, 모든 비분리 사전을 이런 공통 더미로 표현하지는 못합니다.

HMC에서도 방향별 척도를 맞춥니다. 운동량 $p\sim N(0,M)$, 운동에너지 $p^TM^{-1}p/2$라는 규약을 고정하자. 목표가 $N(0,\Sigma)$이면

\[ \dot q=M^{-1}p,\quad \dot p=-\Sigma^{-1}q, \qquad \ddot q=-M^{-1}\Sigma^{-1}q. \]

모든 진동수를 1로 맞추는 선택은 $M=\Sigma^{-1}$입니다. 이 규약에서 $M=\Sigma$라고 쓰면 반대 방향으로 전처리합니다. 소프트웨어가 부르는 inverse mass와 mass를 반드시 구분해야 합니다.

목표 공분산 고윳값 0.01 1 100에서 단위 질량과 정밀도 질량이 만드는 진동수 비교 — 그림 160 목표 공분산 고윳값이 $0.01,1,100$일 때 단위 질량의 진동수는 $10,1,0.1$, $M=\Sigma^{-1}$의 진동수는 모두 1이다. 세로축은 로그 눈금이다. 이 그림은 선형 Hamilton 방정식의 진동수이며 HMC의 유효표본크기 측정은 아니다.#

마지막으로 $\Omega\succ0$의 Gaussian 목적함수 $-\log\det\Omega+\operatorname{tr}(S\Omega)$는 볼록합니다. 대칭 방향 $H$의 두 번째 미분은 $\operatorname{tr}(\Omega^{-1}H\Omega^{-1}H)=\|\Omega^{-1/2}H\Omega^{-1/2}\|_F^2$입니다.

희소 Cholesky가 있으면 $\log\det\Omega=2\sum\log L_{ii}$를 구하되, 일반 희소성은 소거 중 fill-in을 허용합니다. 삼중대각처럼 $O(T)$가 보장되는 구조와 구분합니다. Rademacher 벡터 $z$는 $E[z^TAz]=\operatorname{tr}A$를 주므로 trace 추정에 쓰입니다. $\log\det\Omega=\operatorname{tr}\log\Omega$에 적용할 때는 무작위 추정오차와 $\log\Omega$ 곱의 근사오차를 따로 보고합니다.

X=np.array([[1.,1.]])
V=np.linalg.solve(np.eye(2)+X.T@X,np.eye(2))
Vdual=np.eye(2)-X.T@np.linalg.solve(np.eye(1)+X@X.T,X)
assert np.allclose(V,Vdual)
assert np.allclose(V@X.T@np.array([3.]),[1,1])
# A rectangular factor also handles exact singular covariance.
A=np.array([[2.],[1.]])
assert np.array_equal(A@A.T,[[4,2],[2,1]])
# Complete-pivot Cholesky: residual and numerical rank are explicit.
def pivoted_factor(S,tol=1e-12):
    residual=np.array(S,dtype=float,copy=True)
    cols=[]
    for _ in range(len(S)):
        pivot=int(np.argmax(np.diag(residual)))
        value=residual[pivot,pivot]
        if value<=tol: break
        col=residual[:,pivot]/np.sqrt(value)
        cols.append(col); residual-=np.outer(col,col)
        residual=(residual+residual.T)/2
    if np.linalg.eigvalsh(residual).min() < -tol:
        raise ValueError('input is not PSD within tolerance')
    factor=np.column_stack(cols) if cols else np.zeros((len(S),0))
    return factor,residual
F,residual=pivoted_factor(A@A.T)
assert F.shape==(2,1) and np.allclose(F@F.T,A@A.T)
assert np.linalg.norm(residual,2)<1e-12
# All sign vectors: exact trace identity, not a Monte Carlo claim.
Z=np.array([[1,1],[1,-1],[-1,1],[-1,-1]])
assert np.allclose(np.mean(np.einsum('bi,ij,bj->b',Z,V,Z)),np.trace(V))
print('posterior covariance:',V,'; log precision determinant:',2*np.log(np.diag(L)).sum())
print('pivoted rank:',F.shape[1],'; residual:',np.linalg.norm(residual,2))

posterior covariance: [[ 0.66666667 -0.33333333]
 [-0.33333333  0.66666667]] ; log precision determinant: 2.5649493574615363
pivoted rank: 1 ; residual: 0.0

연습문제와 전체 풀이#

1. 같은 Gaussian을 만드는 두 인수. 첫 $\Sigma$의 Cholesky와 고유분해 인수가 직교변환만큼 다름을 보여라.

풀이. Cholesky는 $L=\left(\begin{smallmatrix}2&0\\1&\sqrt2\end{smallmatrix}\right)$입니다. 고윳값은 $(7\pm\sqrt{17})/2$이며 각 $\lambda$의 고유벡터는 $(2,\lambda-4)$를 정규화하여 얻습니다. 열로 모은 $U$와 $A=U\Lambda^{1/2}$에 대해 $Q=L^{-1}A$라 두면 $QQ^T=L^{-1}\Sigma L^{-T}=I$입니다. 정방행렬이므로 $Q$는 직교이고 $A=LQ$. $Qz\sim N(0,I)$이므로 분포가 같습니다. 특정 $z$에 대한 출력까지 같을 필요는 없습니다.

2. 조건부 독립과 주변 독립. 세 변수 예에서 부분상관 $\rho_{13\cdot2}$와 주변상관 $\rho_{13}$을 구하라.

풀이. 부분상관은 $-\Omega_{13}/\sqrt{\Omega_{11}\Omega_{33}}=0$. 주변상관은 $(1/4)/\sqrt{(3/4)(3/4)}=1/3$. 나머지를 조건으로 한 두 변수 밀도의 교차항은 $-\Omega_{13}x_1x_3$이므로 0일 때 인수분해됩니다. Gaussian 가정 아래에서만 이 계산이 독립 판정입니다.

3. 세 기간 사후의 정규화. 첫 경로 예의 사후 밀도와 둘째 상태의 주변분포를 구하라.

풀이. $p(x\mid y)=\sqrt{13}(2\pi)^{-3/2}\exp[-(x-m)^T\Omega(x-m)/2]$입니다. $x_2\mid y\sim N(5/13,6/13)$. 정밀도 대각 $3$의 역수 $1/3$은 $x_1,x_3$까지 조건으로 주었을 때의 분산이며 주변분산 $6/13$과 다릅니다.

4. 대형 BVAR의 메모리. 절편 없이 $m=40$, 지연 13개인 VAR의 분리 사전과 일반 사전을 비교하라.

풀이. 방정식당 $k=520$, 전체 계수는 $20800$개입니다. 조밀 공분산 원소 수는 $20800^2=432640000$, binary64 저장은 $3461120000$바이트입니다. 두 Kronecker 인수만 저장하면 $520^2+40^2=272000$개, $2176000$바이트입니다. 조밀 Cholesky 비용의 대표항은 $20800^3/3$, 분리 인수는 $(520^3+40^3)/3$입니다. 일반 사전에서도 희소·반복법을 쓰면 조밀 최악비용보다 줄일 수 있으므로 이 비교는 가능한 모든 알고리즘의 하한이 아닙니다. 공통 더미 관측은 분리 구조만 보존합니다.

5. 질량의 규약. $\Sigma=\operatorname{diag}(1/100,100)$에 $M=I,\Sigma,\Sigma^{-1}$를 사용한 진동수를 구하라.

풀이. 진동수 제곱은 $M^{-1}\Sigma^{-1}$의 고윳값입니다. 각각 $(100,1/100)$, $(10000,1/10000)$, $(1,1)$이므로 진동수는 $(10,1/10)$, $(100,1/100)$, $(1,1)$입니다. $M=\Sigma$는 이 규약에서 척도 불균형을 더 키웁니다.

6. trace 추정량의 분산. 대칭 $A$와 독립 Rademacher $z_i$에 대해 $\operatorname{Var}(z^TAz)$를 구하라.

풀이. $z_i^2=1$이므로 $z^TAz=\operatorname{tr}A+2\sum_{i<j}a_{ij}z_iz_j$. 서로 다른 순서 없는 쌍의 곱은 적어도 한 좌표가 홀수 번 나타나므로 기대값이 0입니다. 따라서 분산은 $4\sum_{i<j}a_{ij}^2=2\sum_{i\ne j}a_{ij}^2$. 독립 $s$개 평균의 분산은 이를 $s$로 나눈 값입니다.

지금까지의 내용을 수학의 언어로 정리해 봅시다#

앞에서는 동일한 Gaussian을 공분산, 정밀도, 상태경로라는 세 표현으로 계산했습니다. 이제 표본추출과 조건부 갱신이 왜 정확한지, 어떤 양정치 조건이 필요한지 정리합니다.

정리 1 · Gaussian 인수와 조건부 독립. $z\sim N(0,I_r)$, $AA^T=\Sigma\succeq0$이면 $\mu+Az\sim N(\mu,\Sigma)$입니다. $\Sigma\succ0$이면 $x_i\perp x_j\mid x_{-ij}$와 $\Omega_{ij}=0$이 동치입니다.

증명. 임의의 $t$에 대해 $t^T(\mu+Az)$는 평균 $t^T\mu$, 분산 $\|A^Tt\|^2=t^T\Sigma t$인 정규변수입니다. 이는 특이 경우까지 포함하는 다변량 정규의 정의이며, 특성함수 $\exp(it^T\mu-t^T\Sigma t/2)$로 유일합니다. 조건부 독립은 양정치 밀도에서 확인합니다. 나머지 좌표를 고정하면 두 변수에 관한 로그밀도는 개별 이차항·일차항과 $-\Omega_{ij}x_ix_j$의 합입니다. $\Omega_{ij}=0$이면 두 개의 양의 적분가능 함수의 곱으로 인수분해되므로 독립입니다.

역으로 독립인 양의 매끄러운 밀도의 로그는 각 변수의 함수의 합입니다. 혼합편미분이 0이어야 하므로 $-\Omega_{ij}=0$. $\square$

정리 2 · Gaussian 회귀의 두 갱신식. 위 $V_0,R\succ0$ 가정에서 사후 평균·공분산은 계수공간 식과 관측공간 식 모두로 주어집니다.

증명. 우도와 사전의 지수에 들어가는 제곱을 전개하면 $\beta^TK\beta-2h^T\beta+\text{상수}$입니다. $u^TKu=u^TV_0^{-1}u+(Xu)^TR^{-1}(Xu)>0$이므로 $K$는 양정치입니다. $b=K^{-1}h$를 넣으면 $(\beta-b)^TK(\beta-b)+\text{상수}$여서 평균·공분산이 유일하게 정해집니다.

$S=R+XV_0X^T$라 놓겠습니다. $KV_0X^T=X^TR^{-1}S$이므로

\[ K(V_0-V_0X^TS^{-1}XV_0) =I+X^TR^{-1}XV_0-X^TR^{-1}XV_0=I. \]

이것이 공분산 식입니다. 또 $K[b_0+V_0X^TS^{-1}(y-Xb_0)]=V_0^{-1}b_0+X^TR^{-1}y=h$입니다. $K$가 가역이므로 평균 식도 성립합니다. $\square$

정리 3 · 경로 정밀도와 두 평활 추출법. 독립 Gaussian 초기 상태·전이·관측잡음이 양정치 공분산을 가지면 사후 정밀도는 위 블록 삼중대각이며, 정밀도 추출과 뒤방향 조건부 추출은 동일한 사후분포를 생성합니다.

증명. $\|x_1-m_0\|_{P_0^{-1}}^2+\sum\|x_{t+1}-A_tx_t\|_{Q_t^{-1}}^2+\sum\|y_t-H_tx_t\|_{R_t^{-1}}^2$를 전개합니다. 각 전이항의 두 대각 기여는 $A_t^TQ_t^{-1}A_t,Q_t^{-1}$, 교차 기여는 $-Q_t^{-1}A_t$입니다. 앞서 표시한 블록이 모두 얻어집니다. 동차 이차형식이 0이면 초기항에서 $x_1=0$, 전이항에서 순서대로 $x_2=\cdots=x_T=0$이므로 양정치입니다. 완전제곱으로 $N(\Omega^{-1}h,\Omega^{-1})$를 얻고, 삼각계 추출의 공분산 계산으로 첫 방법의 정확성이 따릅니다.

Markov 밀도는 $p(x_1)\prod p(x_{t+1}\mid x_t)\prod p(y_t\mid x_t)$로 분해됩니다. $x_{t+1}$을 고정하면 이후 관측의 인자는 $x_t$에 의존하지 않습니다. 따라서 $p(x_t\mid x_{t+1:T},y_{1:T})=p(x_t\mid x_{t+1},y_{1:t})$이고, 필터의 결합 Gaussian에 Schur 조건부 공식을 적용하면 본문의 $J_t$ 식입니다. 연쇄법칙으로 $p(x_T\mid y)\prod_{t<T}p(x_t\mid x_{t+1},y_{1:t})=p(x_{1:T}\mid y)$. 두 번째 방법도 같은 정규화된 밀도를 갖습니다. $\square$

삼각 인자와 조건부 분포를 이용하면 같은 확률모형을 여러 경로로 생성할 수 있습니다. 다음 장부터는 유한차원에서 익숙했던 존재와 수렴의 논증을 무한차원 공간에서 다시 점검합니다. I1로 이어 읽기.