11.4 다변량 가우시안

실습 · 연습문제 2개

다변량 가우시안의 조건부분포 condition-gaussian

x∼N(μ,Σ)\mathbf{x} \sim \mathcal{N}(\boldsymbol{\mu}, \boldsymbol{\Sigma})의 성분을 관측하지 않은 것 aa와 관측한 것 bb로 나누면

xa∣xb∼N(μa+ΣabΣbb−1(xb−μb),  Σaa−ΣabΣbb−1Σba)\mathbf{x}_a \mid \mathbf{x}_b \sim \mathcal{N}\big(\boldsymbol{\mu}_a + \boldsymbol{\Sigma}_{ab}\boldsymbol{\Sigma}_{bb}^{-1}(\mathbf{x}_b - \boldsymbol{\mu}_b),\ \ \boldsymbol{\Sigma}_{aa} - \boldsymbol{\Sigma}_{ab}\boldsymbol{\Sigma}_{bb}^{-1}\boldsymbol{\Sigma}_{ba}\big)

이다. 함수 condition(mu, Sigma, obs_idx, x_obs)를 완성한다.

인자 뜻 shape
mu 평균 (D,)
Sigma 공분산 (대칭, 양의 정부호) (D, D)
obs_idx 관측한 성분의 인덱스 목록 길이 k
x_obs 관측값 (obs_idx 순서) (k,)

반환값은 (cond_mean, cond_cov)이고, 관측하지 않은 성분들을 인덱스 오름차순으로 놓은 조건부 평균 (D - k,)와 조건부 공분산 (D - k, D - k)다.

테스트는 다음을 확인한다.

  • 2차원에서 μ2+ρσ2σ1(x1−μ1)\mu_2 + \rho\frac{\sigma_2}{\sigma_1}(x_1 - \mu_1), σ22(1−ρ2)\sigma_2^2(1 - \rho^2)와 같다.
  • 정밀도 행렬 Λ=Σ−1\boldsymbol{\Lambda} = \boldsymbol{\Sigma}^{-1}로 쓴 공식 Σa∣b=Λaa−1\boldsymbol{\Sigma}_{a \mid b} = \boldsymbol{\Lambda}_{aa}^{-1}과 같다.
  • 조건부 공분산은 관측값 x_obs에 의존하지 않는다.
  • 관측 인덱스가 정렬되어 있지 않아도 맞는다.

힌트: Σbb−1\boldsymbol{\Sigma}_{bb}^{-1}를 직접 만들지 말고 np.linalg.solve(Sigma_bb, ...)를 쓴다. 인덱스 배열로 블록을 뽑을 때는 Sigma[np.ix_(a, b)]가 편하다.

선형-가우시안 모델의 사후분포 linear-gaussian-posterior

잠재 상태 x∈RD\mathbf{x} \in \mathbb{R}^D에 대한 사전분포와 선형 관측 모델

x∼N(m0,P0),y=Ax+ε,ε∼N(0,R)\mathbf{x} \sim \mathcal{N}(\mathbf{m}_0, \mathbf{P}_0), \qquad \mathbf{y} = \mathbf{A}\mathbf{x} + \boldsymbol{\varepsilon}, \quad \boldsymbol{\varepsilon} \sim \mathcal{N}(\mathbf{0}, \mathbf{R})

이 있을 때, 사후분포 p(x∣y)=N(m,P)p(\mathbf{x} \mid \mathbf{y}) = \mathcal{N}(\mathbf{m}, \mathbf{P})를 구하는 함수 posterior(m0, P0, A, R, y)를 완성한다.

인자 shape
m0 (D,)
P0 (D, D)
A (M, D)
R (M, M)
y (M,)

반환값은 (m, P)이고 shape은 (D,), (D, D)다.

본문에서 유도한 두 가지 형태 중 어느 것을 써도 된다.

  • gain 형태: S=AP0A⊤+R\mathbf{S} = \mathbf{A}\mathbf{P}_0\mathbf{A}^\top + \mathbf{R}, K=P0A⊤S−1\mathbf{K} = \mathbf{P}_0\mathbf{A}^\top\mathbf{S}^{-1}, m=m0+K(y−Am0)\mathbf{m} = \mathbf{m}_0 + \mathbf{K}(\mathbf{y} - \mathbf{A}\mathbf{m}_0), P=P0−KAP0\mathbf{P} = \mathbf{P}_0 - \mathbf{K}\mathbf{A}\mathbf{P}_0
  • 정밀도 형태: P−1=P0−1+A⊤R−1A\mathbf{P}^{-1} = \mathbf{P}_0^{-1} + \mathbf{A}^\top\mathbf{R}^{-1}\mathbf{A}, m=P(P0−1m0+A⊤R−1y)\mathbf{m} = \mathbf{P}\big(\mathbf{P}_0^{-1}\mathbf{m}_0 + \mathbf{A}^\top\mathbf{R}^{-1}\mathbf{y}\big)

테스트는 다음을 확인한다.

  • 1차원에서 사후 정밀도 = 사전 정밀도 + 측정 정밀도이고, 사후 평균은 정밀도로 가중한 평균이다.
  • 두 형태의 결과가 같다.
  • 관측을 두 번에 나누어(첫 사후를 두 번째의 사전으로) 반영해도, 한 번에 반영해도 같다.
  • 관측이 많아질수록 사후 공분산의 대각 성분이 줄어든다(늘지 않는다).

힌트: 역행렬 대신 np.linalg.solve를 쓰면 수치적으로 더 안정하다. 반환하는 P는 대칭이어야 한다.