메트로폴리스 (Metropolis)

지난 몬테카를로 시뮬레이션 관련 포스트들은 모두 독립 확률 변수에 기초한 것입니다. 독립 확률 변수란 확률 변수들 사이에 상관관계가 없어 각 확률 변수를 독립적으로 만들 수 있다는 말이죠. 수식으로는 두 확률 변수 $X, Y$에 대해 $p(X,Y)=p(X)*p(Y)$와 같이 표현할 수 있을 겁니다. 하지만 확률 변수들을 독립적이지 않을 때는 어떨까요? 다시 말해서 $n$개의 확률 변수 $x_1, ..., x_n$에 대해서 다음 수식처럼 확률 변수를 독립적으로 분해할 수 없다는 것이죠. 이럴 때에는 메트로폴리스 (metropolis) 알고리즘을 사용할 수 있습니다.

$p(x_1, ..., x_n) !=p(x_1)...p(x_n)$

메트로폴리스 알고리즘은 다음과 같이 동작합니다.

1. 정의역에 있는 독립 난수의 집합 $(x_1^1, x_2^1, ..., x_n^1)$ 로부터 시작합니다.
2. 정의역에 있는 다른 독립 난수의 집합 $x^{i+1}=(x_1^{i+1}, ..., x_n^{i+1})$을 만들고 이는 임의의 확률 함수 $Q(x^{i})$를 사용해서 만들 수 있습니다. 함수 $Q$의 유일한 조건은 현재 점 $x$에서 새로운 점 $y$로 갈 확률과 현재 점 $y$에서 새로운 점 $x$로 이동하는 확률이 같아야 한다는 것입니다. (간단하게 균등 함수로 구현할 수 있겠죠)
3. 균등 난수 $z$를 만듭니다.
4. 만일 $p(x^{i+1})/p(x^{i}) < z$ 라면, $x^{i+1}=x^{i}$ 이고 2 단계로 다시 돌아갑니다.
5. 매우 큰 수 $i$에 대해 이 과정을 반복하면 $x^i$의 집합은 목표로 하는 확률 질량 함수 $p(x)$를 가지게 됩니다.

 

파이썬 구현

예를 들어 확률 질량 함수가 $p(x_1, x_2)=6(x_0-x_1)^2$이고 $[0,1]$에 있는 두 난수 $x_1, x_2$를 만들어내는 간단한 경우에 대해 파이썬으로 구현해보도록 하겠습니다.

import random

def metropolis(p, q, x=None):
    while True:
        x_old = x
        x = q(x)
        if x_old is None:
            x_old = [0, 1]
        if p(x) / p(x_old) < random.random():  
            x = x_old
        yield x

def P(x):
    return 6.*(x[0]-x[1])**2

def Q(x):
    return [random.random(), random.random()]

for i, x in enumerate(metropolis(P, Q)):
    print(x)
    if i == 100: break

• x_old 변수는 처음에 x가 None 객체이므로 임의의 값을 지정합니다.
• $Q$는 정의역 $[0,1]\times[0,1]$에서 무작위 점을 만들어내는 함수입니다. return 대신에 생성자 yield 함수를 사용하여 모든 계산을 한 번에 하지 않고 생성된 값으로 반복합니다.

실제로 이 코드를 돌려보면 반복된 값을 많이 만들어냅니다. 즉, 만들어낸 벡터가 이전의 벡터와 두터운 상관관계를 가질 수 있다는 것이죠. 이런 이유로 메트로폴리스 알고리즘에서는 일부를 건너 뛰어서 벡터들간의 상관관계를 없애야 합니다. 이전의 벡터로부터 독립적인 다음 벡터를 만들기 위해 충분히 큰 값을 갖는 decorrelation length 를 설정하고 이 길이에 다다를 때마다 벡터를 샘플링합니다.

def metropolis_decorrelate(p, q, x=None, ds=100):
    k = 0
    for x in metropolis(p, q, x):
        k += 1
        if k % ds == ds -1:
            yield x

for i, x in enumerate(metropolis_decorrelate(P,Q)):
    print(x)
    if i == 100: break