다운로드
작성자: admin 작성일시: 2016-05-25 00:11:51 조회수: 6213 다운로드: 366
카테고리: 기초 수학 태그목록:

다변수 가우시안 정규 분포

$D$차원 다변수 가우시안 정규 분포(MVN: multivariate Gaussian normal distribution)의 확률밀도함수는 평균 벡터 $\mu$ 와 공분산 행렬 $\Sigma$ 라는 두 개의 모수를 가지며 다음과 같은 수식으로 정의된다.

$$ \mathcal{N}(x ; \mu, \Sigma) = \dfrac{1}{(2\pi)^{D/2} |\Sigma| ^{1/2}} \exp \left( -\dfrac{1}{2} (x-\mu)^T \Sigma^{-1} (x-\mu) \right) $$

이 식에서 각 기호의 의미는 다음과 같다.

  • $x \in \mathbf{R}^D $ 확률 변수 벡터
  • $\mu \in \mathbf{R}^D $ 평균 벡터
  • $\Sigma \in \mathbf{R}^{D\times D} $ 공분산 행렬
  • $\Sigma^{-1} = \Lambda \in \mathbf{R}^{D\times D} $ 정밀도 행렬(precision matrix)

이 때 공분산 행렬은 일반적으로 정부호 대칭 행렬(positive definite symmetric matrix)이다.

SciPy의 다변수 정규 분포 명령

SciPy의 stats 서브패키지에는 다변수 정규 분포를 위한 multivariate_normal 명령이 있다. mean 인수로 평균 벡터를, cov 인수로 공분산 행렬을 받는다.

다변수 정규 분포의 예

2차원($D=2$) 다변수 정규 분포의 예를 몇가지 살펴보자.

우선 2차원이므로 확률 변수 벡터는 $$ x = \begin{bmatrix}x_1 \\ x_2 \end{bmatrix} $$ 이다.

경우 1

만약

$$ \mu = \begin{bmatrix}2 \\ 3 \end{bmatrix}. \;\;\; \Sigma = \begin{bmatrix}1 & 0 \\ 0 & 1 \end{bmatrix} $$

이면

$$ | \Sigma| = 1. \;\;\; \Sigma^{-1} = \begin{bmatrix}1 & 0 \\ 0 & 1 \end{bmatrix} $$$$ (x-\mu)^T \Sigma^{-1} (x-\mu) = \begin{bmatrix}x_1 - 2 & x_2 - 3 \end{bmatrix} \begin{bmatrix}1 & 0 \\ 0 & 1 \end{bmatrix} \begin{bmatrix}x_1 - 2 \\ x_2 - 3 \end{bmatrix} = (x_1 - 2)^2 + (x_2 - 3)^2 $$$$ \mathcal{N}(x_1, x_2) = \dfrac{1}{2\pi} \exp \left( -\dfrac{1}{2} \left( (x_1 - 2)^2 + (x_2 - 3)^2 \right) \right) $$

이 확률밀도함수의 모양은 다음과 같다.

In [1]:
mu = [2, 3]
cov = [[1, 0], [0, 1]]
rv = sp.stats.multivariate_normal(mu, cov)
xx = np.linspace(-1, 6, 120)
yy = np.linspace(-1, 6, 150)
XX, YY = np.meshgrid(xx, yy)
plt.contour(XX, YY, rv.pdf(np.dstack([XX, YY])))
plt.axis("equal")
plt.xlim(0, 4)
plt.ylim(2, 4)
plt.xlabel("$x_1$")
plt.ylabel("$x_2$")
plt.title("경우 1의 확률밀도함수")
plt.show()

경우 2

만약

$$ \mu = \begin{bmatrix}2 \\ 3 \end{bmatrix}. \;\;\; \Sigma = \begin{bmatrix}2 & 3 \\ 3 & 7 \end{bmatrix} $$

이면 $$ |\Sigma| = 5,\;\;\; \Sigma^{-1} = \begin{bmatrix}1.4 & -0.6 \\ -0.6 & 0.4 \end{bmatrix} $$

$$ (x-\mu)^T \Sigma^{-1} (x-\mu) = \begin{bmatrix}x_1 - 2 & x_2 - 3 \end{bmatrix} \begin{bmatrix}1.4 & -0.6 \\ -0.6 & 0.4\end{bmatrix} \begin{bmatrix}x_1 - 2 \\ x_2 - 3 \end{bmatrix} = \dfrac{1}{10}\left(14(x_1 - 2)^2 - 12(x_1 - 2)(x_2 - 3) + 4(x_2 - 3)^2\right) $$$$ \mathcal{N}(x_1, x_2) = \dfrac{1}{2\sqrt{5}\pi} \exp \left( -\dfrac{1}{10}\left(7(x_1 - 2)^2 - 6(x_1 - 2)(x_2 - 3) + 2(x_2 - 3)^2\right) \right) $$

이 확률밀도함수의 모양은 다음과 같다.

In [2]:
mu = [2, 3]
cov = [[2, 3], [3, 7]]
rv = sp.stats.multivariate_normal(mu, cov)
xx = np.linspace(-1, 6, 120)
yy = np.linspace(-1, 6, 150)
XX, YY = np.meshgrid(xx, yy)
plt.contour(XX, YY, rv.pdf(np.dstack([XX, YY])))
plt.axis("equal")
plt.xlim(0, 4)
plt.ylim(2, 4)
plt.xlabel("$x_1$")
plt.ylabel("$x_2$")
plt.title("경우 2의 확률밀도함수")
plt.show()

가우시안 정규 분포와 고유값 분해

다변수 가우시안 정규 분포의 공분산행렬 $\Sigma$은 대칭행렬이므로 대각화 가능(diagonalizable)하다. 이 식에서 $\Lambda$는 고유값 행렬, $V$는 고유벡터 행렬이다.

$$ \Sigma^{-1} = V \Lambda^{-1} V^T $$

이를 이용하면 확률밀도함수는 다음처럼 좌표 변환할 수 있다.

$$ \begin{eqnarray} \mathcal{N}(x) &\propto& \exp \left( -\dfrac{1}{2} (x-\mu)^T \Sigma^{-1} (x-\mu) \right) \\ &=& \exp \left( -\dfrac{1}{2} (x-\mu)^T V \Lambda^{-1} V^T (x-\mu) \right) \\ &=& \exp \left( -\dfrac{1}{2} (V^T(x-\mu))^T \Lambda^{-1} (V^T (x-\mu)) \right) \\ &=& \exp \left( -\dfrac{1}{2} (V^{-1}(x-\mu))^T \Lambda^{-1} (V^{-1} (x-\mu)) \right) \\ &=& \exp \left( -\dfrac{1}{2} x'^T \Lambda^{-1} x' \right) \\ \end{eqnarray} $$

즉, 고유벡터를 기저벡터로 사용하여 $x$를 $x' = V^{-1}(x-\mu)$ 로 좌표 변환하면, 좌표 변환된 새로운 확률변수 $x'$의 공분산행렬은 대각 행렬인 고윳값 행렬 $\Lambda$가 된다. 공분산 행렬이 대각행렬이므로 좌표 변환된 $x'$은 서로 독립이다.

예를 들어 위의 두번째 예제에서 공분산행렬을 고유 분해하면 다음과 같다.

In [3]:
mu = [2, 3]
cov = [[4, 3], [3, 5]]
w, V = np.linalg.eig(cov)
In [4]:
w
Out:
array([1.45861873, 7.54138127])
In [5]:
V
Out:
array([[-0.76301998, -0.6463749 ],
       [ 0.6463749 , -0.76301998]])

확률밀도함수가 고유벡터 $v_1 = (-0.763, 0.646)$와 $v_2 = (-0.646, -0.763)$를 축으로하는 타원형임을 알 수 있다.

In [6]:
xx = np.linspace(-1, 5, 120)
yy = np.linspace(0, 6, 150)
XX, YY = np.meshgrid(xx, yy)

plt.figure(figsize=(8, 4))

d = dict(facecolor="k", edgecolor="k")

plt.subplot(121)
rv1 = sp.stats.multivariate_normal(mu, cov)
plt.contour(XX, YY, rv1.pdf(np.dstack([XX, YY])))
plt.annotate("", xy=(mu + 0.35 * w[0] * V[:, 0]), xytext=mu, arrowprops=d)
plt.annotate("", xy=(mu + 0.35 * w[1] * V[:, 1]), xytext=mu, arrowprops=d)
plt.scatter(mu[0], mu[1], s=10, c="k")
plt.axis("equal")
plt.xlabel("$x_1$")
plt.ylabel("$x_2$")
plt.title("$X_1$,$X_2$의 결합 확률밀도함수")

plt.subplot(122)
rv2 = sp.stats.multivariate_normal(mu, w)  # 좌표 변환
plt.contour(XX, YY, rv2.pdf(np.dstack([XX, YY])))
plt.annotate("", xy=(mu + 0.35 * w[0] *
                     np.array([1, 0])), xytext=mu, arrowprops=d)
plt.annotate("", xy=(mu + 0.35 * w[1] *
                     np.array([0, 1])), xytext=mu, arrowprops=d)
plt.scatter(mu[0], mu[1], s=10, c="k")
plt.axis("equal")
plt.xlabel("$x'_1$")
plt.ylabel("$x'_2$")
plt.title("$X'_1$,$X'_2$의 결합 확률밀도함수")

plt.tight_layout()
plt.show()