그림으로 깁스샘플링 이해하기

깁스 샘플링에 대해서는 토픽 모델링 공부를 하면서 귀에 피가 나도록 들었었는데요, 사실 그 실체를 정확하게 파악하는데에는 꽤 오랜 시간이 걸렸었습니다. 막연하게만 이해하고 있던 깁스 샘플링을 정확하게 깨닫게 된 건 어떤 2장 짜리 논문에서 본 그림 덕분이었는데, 이번 포스팅에서는 그걸 공유해보고자 합니다.

참고한 논문은 다음과 같습니다.

Breslaw, J. A. (1994). Random sampling from a truncated multivariate normal distribution. Applied Mathematics Letters, 7(1), 1-6.

다변수 정규 분포

다변수 정규 분포(Multivariate Normal Distribution)에 대해서는 상관 토픽 모델 포스팅에서 살짝 설명했었는데요, 우리가 흔히 보던 종 모양의 정규분포 곡선을 2차원 이상으로 확장한 것이라고 보면 됩니다. 변수가 2개 이상이 되면, 단일 변수의 분산뿐만 아니라 변수들 간의 분산(공분산)이라는 개념도 생겨나게 됩니다. 이는 두 변수가 얼마나 같이 움직이는지를 보여주는 척도가 됩니다. 아래의 그래프는 2차원 정규분포의 밀도를 표현한 것인데, X와 Y의 평균은 둘 다 0, X와 Y의 분산도 둘다 1이고, X와 Y의 공분산은 0.8인 경우입니다. 괜시리 어렵게 수식으로 쓰면 다음과 같겠네요.

위 그림은 [-2, 2] 범위를 0.1 단위로 격자를 나눈 뒤, 앞서 말한 2차원 정규분포에서 점을 1만개를 뽑아 어느 격자에 속하는지를 센 뒤 그 비율에 따라 색깔을 달리하여 칠한 것입니다. 당연히 중심인 0,0 부분이 제일 진하고 주변으로 갈수록 옅어지는 모양이 됩니다. 사실 [-2, 2] 범위 바깥에서 뽑힌 점들도 있지만, 그 수가 너무 적기 때문에 그림에서는 표현하지 않았습니다.

우리가 이 포스팅에서 관심 가지는 것은 '저런 분포가 주어졌을때 어떻게 거기서 표본을 뽑아낼 것인가'입니다. 다행히도 정규분포에서 임의의 표본을 뽑아내는 방법은 박스 뮬러 변환(Box-Muller Transform)이라는 효과적인 기법이 개발되어 있어서, 이를 이용해 누구나 정규분포에서 손쉽게 표본들을 샘플링할 수 있습니다. 다변수 정규분포 역시 본질적으로 정규분포를 여러 차원에 대해 반복하는 것으로 샘플링할 수 있습니다.

이제 한 발 더 나아가봅시다. 위의 빨간 네모는 -2<X<-1, 0<Y<1인 영역을 표현하는데요, 정규분포 중 저 안에 들어가는 녀석들에 대해서만 표본을 뽑는 문제를 생각해볼 수 있겠습니다.

절단 정규 분포와 기각 샘플링

절단 정규 분포(Truncated Normal Distribution)는 정규 분포에서 특정영역만 잘라내어 만든 분포를 가리킵니다. 다르게 말하면 원래는 X의 범위가 무한했던 정규분포에서 X가 가질 수 있는 범위를 제한시킨 분포라고 할 수 있겠네요. 종종 복잡한 확률 과정을 계산하다보면 정규 분포의 정의역을 여러 가지로 제한하게 되는데 이때 사용하게 되는 분포입니다. 이번 예제에서는 위의 2변수 정규분포에서 정의역을 -2<X<-1, 0<Y<1로 제한(빨간색 박스)한 절단 정규 분포를 생각해보겠습니다. 

위의 빨간 박스부분만 떼어서 다시 그림을 그리면 아래와 같아집니다. 원래 -2<X<-1, 0<Y<1 범위 내의 표본 수는 매우 적었지만, 이제 전체 구간을 줄여버렸으므로 상대적으로 (-1, 0) 부분의 색이 진해집니다. 이제 이 범위 내에 속하는 표본을 뽑는 작업을 생각해봅시다. 이는 생각보다 어려운 작업인데요, 왜냐하면 전체 2차원 정규분포에서 표본을 뽑으면 이 구간에 들어가지 않을 확률이 들어갈 확률보다 훨씬 높기 때문입니다. (위의 2차원 정규분포의 중심은 (0,0)에 있었다는 걸 기억해보세요. 대부분의 표본은 그 근처에서 뽑힐 겁니다.)

그래도 끈기를 가지고 수많은 표본을 뽑다보면 우연히 우리가 원하는 사각형 내에 들어오는 샘플이 있을테니깐 걔네를 열심히 모아보도록 합시다. 이 방법을 기각 샘플링(Rejection Sampling)이라고 합니다. 이 방법의 장점은 '별 다른 고민없이 기존에 사용하던 샘플링 알고리즘을 그대로 사용할 수 있다'는 것이고, 또한 '의도한 분포대로 표본을 뽑을 수가 있다는 것'입니다. 다만 단점은 '표본의 기각률이 높을 경우 실용적으로 써먹기엔 틀렸다'는 것입니다. 

지금 예시로 든 이 절단 정규분포에서 1만개의 표본을 뽑기 위해서는 원 2변수 정규분포에서 약 188만개의 표본을 뽑아야 합니다. 그래야 그중 1만개 정도가 우리가 원하는 영역에 들어온다는 것이죠. 만약에 2차원이 아니라 3차원 절단 정규분포라면? 기각률은 훨씬 더 높아질 겁니다. 보통 우리가 다루는 분포는 최소 몇십 차원에서 몇 백 차원을 넘나든다는 것을 떠올려보면, 기각 샘플링은 절단 정규분포를 표집하는데에는 적절하지 않다는 결론을 얻게 됩니다.

-보너스-

기각 샘플링을 파이썬 코드로 작성하면 다음과 같습니다.

import numpy as np import matplotlib.pyplot as plt R = 0.8 NUM = 10000 np.random.seed(777) mat = np.zeros((10, 10)) # 넉넉하게 2백만개 뽑기로 했습니다. points = np.random.multivariate_normal([0, 0], [[1, R], [R, 1]], 2000000) points = np.array([(x, y) for x, y in points if -2 <= x < -1 and 0 <= y < 1]) # 조건을 통과한 것중 1만개만 선택합니다. points = points[:NUM] for x, y in points: mat[int(y * 10), int((x + 2) * 10)] += 1 mat /= mat.sum() fig, ax = plt.subplots() ax.matshow(mat, cmap=plt.cm.Blues) plt.show()

각개격파!

그래서 다른 방법을 찾아야 합니다. 우리가 뽑을 표본은 2차원 상의 점이므로, 이를 하나씩 분리해서 처리하면 조금 더 쉽지 않을까요? 사실 앞에서 밝히지 않은 게 있는데, 다변수 정규 분포를 뽑는 과정이 정규 분포를 하나씩 뽑아서 모으는 과정이라는 겁니다. 이걸 그대로 절단 정규 분포에 적용하면 되지 않을까요?

두 변수의 공분산을 ρ이라 할때 다음 2변수 정규분포는 다음과 같이 정규분포 두 개로 분해할 수 있습니다.

일때

라고 하면

가 된다.

여기서 Z1, Z2는 표준 정규분포에서 뽑은 표본입니다. X와 Y는 Z1과 Z2의 조합으로 표현할 수 있습니다. 여기서는 먼저 Z1을 X로 잡았고, 이 Z1의 값에 따라 Y값이 영향을 받게 됩니다. 혹은 정반대로

로 정의해도 무방합니다.

이제 여기서 절단 조건을 추가해보도록 합시다. a<X<b여야 하고, c<Y<d여야 한다고 해봅시다. 그럼 먼저

,

따라서 Z1은 truncN(a, b)에서 표집하면 되고, Z2는 Z1를 표집한 다음 그 값을 이용해 truncN(~~, ~~)에서 표집하면 되겠습니다. (수식 적기 귀찮아서 ~~로 대체했습니다. 다들 생략된게 뭔지 아시죠?) 여기서 truncN은 1변수 절단 정규 분포를 가리킵니다. 다행히도 1변수 절단 정규분포의 경우 샘플링 기법이 발명되었기 때문에 손쉽게 샘플링을 실시할 수 있습니다. 만약 이를 사용할 수 없는 경우라면 정규분포에서 기각 샘플링을 해도 좋습니다. (1변수 정규분포의 경우 다변수에서처럼 기각률이 높지 않기 때문에 사용해봄직합니다.)

파이썬 코드로 구현해보면 다음과 같습니다.

import numpy as np import matplotlib.pyplot as plt R = 0.8 NUM = 10000 np.random.seed(777) def meth1(): x = truncnorm.rvs(-2, -1) std = (1 - R ** 2) ** 0.5 y = R * x + truncnorm.rvs((0 - R * x) / std, (1 - R * x) / std) * std return x, y mat2 = np.zeros((10, 10)) points = np.array([meth1() for _ in range(NUM)]) for x, y in points: mat2[int(y * 10), int((x + 2) * 10)] += 1 mat2 /= mat2.sum() fig, ax = plt.subplots() ax.matshow(mat2, cmap=plt.cm.Blues) plt.show()

자, 이렇게해서 샘플링을 하면 어떤 분포가 나올까요? 짜잔~

모양이 좀 다르죠? 원래 우리가 원했던 모양은 오른쪽 아래만 새파랗고 나머지는 하얘야하는데, 이 방법으로 추출한 경우 Y = 0쪽에 뭔가 전반적으로 분포가 쏠려있습니다. 왜 이렇게 잘못된 분포가 나오게 되었을까요? 이는 우리가 X를 먼저뽑고 그에 따라 조건부 확률로 Y를 뽑았기 때문입니다. 반대로 Y를 먼저 뽑고 X를 Y에 대한 조건부로 뽑았다면 X = -1 쪽에 쏠린 분포가 나오게 됩니다.

어디가 잘못되는지 감이 안온다면 문제를 조금 단순화시켜봅시다. X, Y ∈ {a, b, c}인 확률 변수이고,  a, b, c가 뽑힐 확률은 각각 1:2:3의 비율을 이룬다고 가정해봅시다. 그런데 X와 Y를 뽑는 사건은 독립이 아니라서 X와 Y를 동시에 뽑는 경우는 단순히 X를 뽑을 확률과 Y를 뽑을 확률을 곱해서 구할 수 없고 다음과 같이 계산된다고 가정해봅시다.

	X=a	X=b	X=c
Y=a			0.1
Y=b		0.1	0.2
Y=c	0.1	0.2	0.3

자 만약 이 분포를 근사하기 위해 X를 먼저 뽑고, 그 다음 거기서 Y를 뽑는다면 어떻게 될까요? 먼저 X가 뽑힐 확률은 다음과 같이 됩니다.

	X=a	X=b	X=c
	0.167	0.333	0.5

보면 벌써부터 틀렸습니다. X와 Y의 결합분포에서는 X=a가 되는 경우는 0.1밖에 없으나 Y를 무시하고 X를 뽑으면 X=a가 뽑힐 확률이 0.167로 0.1보다 훨씬 큽니다. X=b인 경우도 마찬가지이죠. X=a와 X=b가 원 분포보다 더 많이 뽑히니 이 상태에서 Y가 뽑힐 확률을 고려해 X와 Y의 결합분포를 계산해보았자 틀린 결과가 나오게 됩니다.

이게 바로 위 그림에서 발생한 현상과 동일한 현상입니다. 원래는 X=-2 근방에서 뽑힐 확률이 더 적어야하지만, 우리는 Y를 무시하고 X에 대해서 우선 샘플링한 뒤 Y를 샘플링하게 되므로, X=-2 근방에서 뽑힐 확률이 증가하여 아래쪽이 두터운, 잘못된 모양의 분포를 산출하게 된 것이지요.

따라서 각개격파 방법은 다변수 절단 정규분포를 샘플링하는데에는 쓸수 없다는 사실을 깨닫게 됩니다. 그렇다면 우리가 원하는 분포를 추출하기 위해서는 기각 샘플링 삽질을 할 수 밖에 없는걸까요?

깁스 샘플링!!!

이럴 때 사용할 수 있는게 깁스 샘플링(Gibbs Sampling, 기브스 표집)입니다. 깁스 샘플링의 아이디어를 한 줄로 줄여 말하면, X와 Y가 서로 종속적인 관계에 있으니깐, X값을 정한상태에서 Y를 뽑아서 Y값을 구하고, 다시 거꾸로 이 Y값에서 X를 뽑아서 X값을 정하자는 것이죠. 수식으로 멋지게 쓰면 다음과 같겠습니다.

풀어보자면, 첫번째 X는 그냥 독립분포인 P(X)에서 뽑고, Y는 뽑은 X값을 조건부로하는 Y분포에서 뽑습니다. 그리고 각자 이전의 X 혹은 Y값을 바탕으로 조건부확률을 계산하여 다음 X혹은 Y값을 뽑는 것이죠. 포인트는 현재 샘플링하는 변수 이외의 나머지 변수들이 조건으로 고정되어 있다고 가정한다는 것입니다. (사실 초기값인 X0, Y0은 꼭 P(X)나 P(Y|X=X0)에서 뽑지 않고 아무 랜덤값을 주어도 됩니다.) 이렇게 새로운 X와 Y를 뽑는 과정을 계속 반복하다보면 여기서 나오는 X_i, Y_i값들이 원래 우리가 구하고자했던 X, Y의 결합 분포에서 추출한 표본이 된다는 겁니다. 신기하죠?

깁스샘플링은 변수가 더 많을 때에도 마찬가지로 잘 작동하기 때문에 변수가 굉장히 많고 그 조건들이 까다로워 전체의 결합 확률분포를 계산해 거기서 샘플을 추출하기 어려울 때 아주 유용하게 사용할 수 있습니다. 이번 포스팅에서 예시로 들고 있는 다변수 절단 정규 분포가 딱 그에 맞는 상황입니다. 조건부 확률에서는 쉽게 샘플링할 수 있지만, 결합 확률에서는 샘플링이 어려운 경우 말입니다.

이 과정을 파이썬 코드로 정리하면 다음과 같습니다.

import numpy as np import matplotlib.pyplot as plt R = 0.8 NUM = 10000 np.random.seed(777) def meth2(sampling=5): x = truncnorm.rvs(-2, -1) std = (1 - R ** 2) ** 0.5 y = R * x + truncnorm.rvs((0 - R * x) / std, (1 - R * x) / std) * std # 샘플링을 여러번 반복합니다! for _ in range(sampling): x = R * y + truncnorm.rvs((-2 - R * y) / std, (-1 - R * y) / std) * std y = R * x + truncnorm.rvs((0 - R * x) / std, (1 - R * x) / std) * std return x, y mat2 = np.zeros((10, 10)) # 일단 테스트로 깁스 샘플링을 1회만 돌려봅시다. points = np.array([meth2(1) for _ in range(NUM)]) for x, y in points: mat2[int(y * 10), int((x + 2) * 10)] += 1 mat2 /= mat2.sum() fig, ax = plt.subplots() ax.matshow(mat2, cmap=plt.cm.Blues) plt.show()

짜잔! 그러면 이렇게 처음에 우리가 그렸던 그림과 거의 동일한 분포를 얻을 수 있습니다. 이 그림은 깁스 샘플링을 1회 돌려서 나온 X, Y값을 사용한 경우이구요, 더 많은 횟수를 돌린 경우에도 위와 거의 동일한 모양의 분포가 나옵니다. 보통 깁스 샘플링 사용시 초기 몇 십 회 이내의 표본은 편향되어 있을 가능성이 높아서 사용하지 않습니다만 (이렇게 앞부분의 샘플들을 생성하는 과정을 Burn-in iteration이라고 합니다. 본격적으로 샘플을 뽑기에 앞서 불을 때는 과정이라고 생각하면 되겠죠), 예시로 든 2변수 절단 분포의 경우 차원수가 낮기 때문에 1회만으로도 원 분포와 거의 유사한 분포가 나왔습니다.

그리고 위 코드에서는 매번 깁스 샘플링 결과를 얻기 위해 초기값으로 X, Y를 설정하고 burn-in을 거쳐 최종 샘플을 뽑는 과정을 반복하고 있는데, 굳이 그럴 필요 없이 매 iteration 마다 생성되는 표본들을 연속적으로 사용해도 동일한 분포를 얻을 수 있습니다. 단, 이 경우 이전의 표본과 바로 다음의 표본은 서로 직접적으로 종속관계에 있으므로 대게 1~2개씩 건너뛰면서 사용하는 경우가 많습니다. (10회째, 12회째, 14회째 ... 결과를 사용하는 것처럼)

수치로 비교해보기

그림으로만 보면 우리가 근사한 분포가 얼마나 원 분포와 얼마나 가까운지 확인하기 어렵습니다. 따라서 두 확률분포의 거리를 측정하는데 쓰이는 KL Divergence 척도를 이용해 유사도를 확인해보았습니다.

	각개격파	깁스 샘플링 1회	3회	5회
KL Div	0.2149	0.0139	0.0130	0.0105

KL Divergence는 두 분포가 완전히 같은 경우 0, 다를수록 큰 양수의 값을 가집니다. 처음 시도했던 각개격파 방법의 경우 0.2149로 매우 큰 오차를 보였지만, 깁스 샘플링의 경우 0.01 정도로 거의 원 분포와 동일한 근사 분포를 보이고 있음을 확인할 수 있습니다. 특히 그 차이는 샘플링 횟수가 늘어날수록 줄어들고 있구요.

정리해보자면, 변수가 많고 복잡한 확률분포에서 표본을 표집하는 것은 때때로 어려운데, 이를 단순한 확률 분포 상의 표집의 반복으로 바꿔주는 장치가 깁스 샘플링이라 할 수 있겠습니다. 본 포스팅에서는 1변수 절단 정규분포를 이용해서 다변수 절단 정규분포를 표집하는 것을 실험해보았구요. 현업에서 사용되는 다양한 베이지안 모델은 이보다 더 복잡한 경우가 많으니, 당연히 이 곳에서도 매우 유용하게 쓰이고 있습니다. 이번 기회에 그 개념을 확실히 잡아서 다음번에 깁스 샘플링을 만나도 당황하지 않으시면 좋겠네요!