무작위 걸음(랜덤 워크 random walk) 103 : 주방정식(Master equation)을 이용한 시간에 따른 랜덤 워크의 확률 밀도 함수 계산

지난 포스팅에서 랜덤 워크를 정의하고, 앙상블 평균의 개념을 도입하여 특정한 값을 정의하고 계산하는 방법에 대해서 알아 보았습니다. 또한 확률 밀도 함수를 정의하고, 확률 밀도 함수로 부터 경로의 함수로 주어지는 특정한 함수의 앙상블 평균을 구하는 방법, 또한 확률 밀도 함수의 스텝에 따른 변화를 기술해 주는 주방정식에 대해서도 설명하였습니다. 지난 포스팅을 읽지 않으신 분이라면, 지난 포스팅을 먼저 읽고 이 다음 내용을 읽으시길 권장합니다. 

https://studyingrabbit.tistory.com/89

무작위 걸음(랜덤 워크 random walk) 101 : 파이썬 프로그램으로 구현

https://studyingrabbit.tistory.com/90

무작위 걸음(랜덤 워크 random walk) 102 : 앙상블 평균과 몬테 카를로 시뮬레이션

이번 포스팅에서는 주방정식을 활용하여 확률 밀도 함수의 변화에 대해서 좀 더 알아 보도록 하겠습니다. 확률 밀도 함수는 그 자체로도 매우 중요할 뿐 아니라, 랜덤 워크 스텝(의 증가)에 따른 확률 밀도 함수의 변화에 대한 직관을 구체적인 경우의 예시를 통해서 키우는 것은 랜덤 워크 혹은 더 나아가 확률 과정(stochastic process)에 대해 공부하는데 필수적인 요소이기 때문입니다. 

 

확률 밀도 함수와 주방정식 복습

 

확률 밀도 함수는 특정한 위치, 특정한 스텝에서 입자(랜덤 워커)를 발견할 확률로 정의 됩니다. 혹은 특정한 위치, 특정한 스텝에 입자(랜덤 워커)가 존재할 확률이라고도 볼 수 있습니다. 여기서의 확률은 앙상블 평균으로 정의되는데,

$$p(x, n) = \sum_{i} \delta_{x, X^i_N(n)} P(X^i_N) = \langle \delta_{x, X(n)} \rangle$$

와 같습니다. 모든 가능한 랜덤 워크 경로 중에서 $n$ 번째 스텝에서 $x$에 위치하는 경로의 비율을 뜻 합니다. $p(x, n+1)$과 $p(x', n)$의 관계에 대한 식을 주방정식이라고 합니다. 대부분의 물리 법칙은 보통 시간 $t$에서 시스템의 상태와 시간 $t + \Delta t$에서 시스템의 상대에 대한 관계를 표현하는데(따라서 대부분의 물리 법칙은 미분방정식으로 주어집니다), 주방정식 역시 $n$번째 스텝에서 시스템의 상태와 $n+1$번째 스텝에서 시스템의 상태에 대한 관계를 나타냅니다. 1차원 격자 공간 랜덤 워크의 경우 주방정식은

$$p(x, n+1) = \frac{1}{2}p(x-1, n) + \frac{1}{2} p(x+1, n)$$

으로 주어집니다. 초기 상태($n=0$)에서의 확률 밀도 함수 $p(x, 0)$이 주어져 있다면, 주방정식을 통해서 임의의 스텝에서 임의의 위치의 확률 밀도 함수를 구할 수 있습니다. 

 

주방정식의 선형성

 

서로 다른 두 확률 밀도 함수 $p_1(x, n)$과 $p_2(x,n)$가 각각 주방정식을 만족한다면, 이 둘을 합친 확률 밀도 함수 역시 주방정식을 만족합니다. 여기서두 확률 밀도 함수의 합은 역시 확률 밀도 함수가 되어야 하기 때문에, 두 확률 밀도 함수의 합을

$$p(x, n) = \alpha p_1(x, n) + (1-\alpha) p_2(x,n)$$

으로 정의할 수 있습니다. 여기서 $0 \le \alpha \le 1$ 입니다. $p(x,n)$이 주방정식의 해가 됨을 보이는 것은 어렵지 않은데,

$$\alpha p_1(x, n+1) = \alpha ( \frac{1}{2}p_1(x-1, n) + \frac{1}{2}p_1(x+1, n)) $$

$$(1-\alpha) p_2(x, n+1) = (1-\alpha) ( \frac{1}{2}p_2(x-1, n) + \frac{1}{2}p_2(x+1, n)) $$

를 더하면, 좌변은 $p(x,n+1)$이 되고, 우변에서는 $\frac{1}{2}$를 앞으로 빼고, $x-1$과 $x+1$가 관여된 항 끼리 각각 더하면,

$$\frac{1}{2} (\alpha p_1(x-1, n) + (1-\alpha)p_2(x-1,n)) + \frac{1}{2} (\alpha p_1(x+1, n) + (1-\alpha)p_2(x+1,n))$$

가 됩니다. 이 값은 정의에 따라서 $\frac{1}{2} p(x-1,n) + \frac{1}{2}p(x+1, n)$과 같고, 따라서 $p(x, n)$ 역시 주방정식의 해가 됩니다. 

 

방정식이 해에 대해서 선형성을 갖고 있다면, 해를 구하거나 방정식의 해를 유추할 때 매우 편리 합니다. 복잡한 초기 조건이 주어졌을 때, 만일 이 초기 조건을 여러개의 쉬운 초기 조건의 합으로 표현할 수 있다면, 각 쉬운 초기 조건에 대해서 해를 먼저 구하고, 각각 구한 해를 전체 더하면 이 해는 복잡한 초기 조건에 대한 해가 되기 때문입니다. 

 

구체적인 초기 확률 밀도 함수에 대해서 주방정식을 적용하여 임의의 스텝에서 확률 밀도 함수를 구해 보도록 하겠습니다. 

 

벡터와 행렬을 이용한 표현

 

주방정식은 $p(x,n+1)$은 $p(x',n)$의 선형합으로 표현되는데, 따라서 이를 벡터와 행렬을 이용하면 매우 깔끔하게 표현할 수 있습니다. $x$는 정수값을 갖기 때문에, $p(x,n)$에 대응되는 벡터 $\vec{p}_n$를

$$\vec{p}_n(x) = p(x,n)$$

와 같이 정의 하도록 하겠습니다. $x$는 양수 뿐 아니라, 0이나 음수값을 가질 수 있는데, 

$$\vec{p}_n = (... p_{-2, n}, p_{-1, n}, p_{0, n}, p_{1, n}, p_{2, n}, ...)$$

로 표현하도록 하겠습니다. "양쪽 방향으로 뻗어나가는 무한차원 벡터"를 생각하면 됩니다. 이에 대한 수학적인 엄밀함은 차치해 두고, 직관적인 방식으로 이를 이해하도록 하겠습니다. 이 벡터 표현법을 통해 $p(x,n)$을 표현하면, 주방정식은

$$\vec{p}_{n+1} =A \vec{p}_{n}$$

여기서, 행렬 $A$는 

$$\begin{equation}
      A_{a,b}=\begin{cases}
        \frac{1}{2},  & |a-b|=1 \\
        0, & \text{otherwise}
      \end{cases}
    \end{equation}$$

가 됩니다. 행렬 $A$의 일부분을 직접 써 보면, 

$$A =  \begin{pmatrix}
\cdot \cdot \cdot &\cdot \cdot \cdot&\cdot \cdot \cdot&\cdot \cdot \cdot&\cdot \cdot \cdot\\
\cdot \cdot \cdot&0&\frac{1}{2}&0&\cdot \cdot \cdot\\
\cdot \cdot \cdot&\frac{1}{2}&0&\frac{1}{2}&\cdot \cdot \cdot\\
\cdot \cdot \cdot&0&\frac{1}{2}&0&\cdot \cdot \cdot\\
\cdot \cdot \cdot &\cdot \cdot \cdot&\cdot \cdot \cdot&\cdot \cdot \cdot&\cdot \cdot \cdot
\end{pmatrix} $$

와 같습니다. 행렬 원소의 대부분은 0인데, 대각성분에서 행이나 열 방향으로 1만큼 차이나는 성분에서의 원소의 값은 $\frac{1}{2}$ 이 됩니다. 행렬 $A$를 이와 같이 정의한다면, 주방정식은 $\vec{p}_{n+1} =A \vec{p}_{n}$로 표현될 수 있음을 간단하게 확인할 수 있습니다. 

 

$n=0$에서 확률 밀도 함수에 대응되는 벡터 $\vec{p}_0$이 주어진다면, $n$번째 스텝의 확률 밀도 함수에 대응되는 벡터 $\vec{p}_n$은 행렬 $A$의 거듭제곱으로 나타낼 수 있는데, 

$$\vec{p}_n = A^n \vec{p}_0$$

이 됩니다. 

 

 

1. $p(0, 0) = 1$

 

가장 간단한 초기 확률 밀도 함수 입니다. 원점에 모든 확률이 존재합니다. 이는 지금까지 우리가 해 왔던, 원점에서 시작하는 랜덤 워크와 동일 합니다. 주방정식에 따라서 임의의 스텝 $n$에서 $p(x,n)$를 구하면 이항 분포의 성질에 따라서,

$$\begin{equation}
      p(x,n)=\begin{cases}
        \frac{1}{2^n} B(n, \frac{n+x}{2}),  & \text{if } \frac{n+x}{2} \text{ is an integer} \\
        0, & \text{otherwise}
      \end{cases}
    \end{equation}$$

가 됩니다. 여기서 $B(a,b)$는 이항 계수 입니다. 

$n$에 따른 $p(x,n)$를 그래프로 표현하면 위와 같습니다. $p(x,9)$의 경우, $x=9, x=-9$에서의 $p(x,9)$값이 $0$인 것 처럼 보이는데, 실제로 그 값은 $0$이 아닌 $\frac{1}{2^9}$ 입니다. 이미 첫 포스팅에서 부터 많이 다룬 것 처럼 스텝 $n$이 커질 수록 원점에서 점점 퍼져나가는 개형을 보여줍니다. 만일 $n$이 크다면, 위 분포는 정규 분포를 닮아가고 정규 분포의 표준편차는 $\sqrt{n}$ 이 됩니다. 

 

2. $p(-3,0) = p(3,0) = \frac{1}{2}$

 

이번에는 $x=3, x=-3$에 절반씩의 확률 밀도 함수가 정의된 초기 상태에서의 $p(x,n)$을 계산해 보도록 하겠습니다. 이 값은 주방정식의 행렬표현 $\vec{p}_n = A^n \vec{p}_0$을 이용하여 바로 계산해도 되고, 아니면 주방정식의 선형성을 이용하여 $p(-3,0) = 1$인 경우의 해 $p(x,n)_{-3}$과 $p(3,0) = 1$인 경우의 해 $p(x,n)_3$을 각각 구하고 이 둘을 더하면 됩니다. $p(-3,0) = 1$의 초기 조건의 해는 $p(0,0)=1$의 초기 조건의 해에서 $x$를 $3$만큼 평행이동 하면 됩니다. 즉, $p(0,0)=1$의 초기 조건의 해를 $p(x,n)_0$이라고 한다면, $p(x,n)_{-3} = p(x-3,n)_0$이 됩니다. 

$p(x,n)$을 그래프로 그리면 위와 같습니다. $x=3, -3$에서 각각 확률 밀도 함수가 펴저 나가는 것의 선형 합이 됩니다. 

 

3. 2차원 격자 공간에서의 주방정식

 

이번에는 차원을 하나 높이겠습니다. 2차원을 다룰 텐데, $n$번째 스텝에서의 2차원 격자점 $(x,y)$에서 확률 밀도 함수를 $p(x,y,n)$이라고 하면, 주방정식은

$$p(x,y,n+1) = \frac{1}{4} \Big( p(x+1, y, n) + p(x-1, y, n) + p(x, y+1, n) + p(x, y-1, n)\Big)$$

로 주어집니다. $(x,y)$ 도달하기 위해서는 바로 전 스텝에서 $(x,y)$의 상/하/좌/우에 있어야 하고, 이동할 수 있는 방향이 네 방향이기 때문에 평균적으로 $\frac{1}{4}$에 해당하는 확률 밀도가 각 방향으로 이동하기 때문에 괄호 앞에 $\frac{1}{4}$가 오게 됩니다. 벡터와 행렬을 이용한 표기법 역시 1차원인 경우와 동일한데, 행렬 $A$를 편의상 $A(\vec{x}, \vec{x}')$이라고 쓴다면, 

$$\begin{equation}
      A(\vec{x}, \vec{x}')=\begin{cases}
        \frac{1}{4},  & |\vec{x} - \vec{x}'| = 1 \\
        0, & \text{otherwise}
      \end{cases}
    \end{equation}$$

가 됩니다. 여기서 $\vec{x}$는 당연히 $(x,y)$를 의미합니다. 

초기 조건 $p(0,0,0) = 1$ 에 따른 확률 밀도 함수의 변화는 위와 같습니다. 위 애니메이션에서 원의 넓이는 $p(x,y,n)$값의 크기에 비례합니다. 원점에 국한돼 있던 확률 밀도 함수는 스텝의 증가에 따라 주변부로 퍼지게 됩니다. 

 

임의의 $n$에 대해서 $p(x, y, n)$을 알고 있기 때문에 이 값을 이용해서 임의의 $x, y$의 함수에 대한 기대값(앙상블 평균)을 구할 수 있습니다. 함수 $F(x,y)$가 주어졌을 때, $n$ 번째 스텝에서 $F$의 랜덤 워크의 앙상블 평균은

$$\langle F \rangle(n) = \sum_{x, y} F(x, y) p(x, y, n)$$

입니다. 이 식을 이용하여 $\sigma = \sqrt{\langle r^2 \rangle}(n)$을 구하면,

(위 그래프의 $y$축은 $\sqrt{\langle r^2 \rangle}$ 입니다. )

 

결과는 위 그래프와 같습니다. 랜덤 워크의 스텝 $n$은 정수이기 때문에, 실제로 구할 수 있는 값은 푸른색 점으로 표시한 값이고, 추세를 표현하기 위해서 $\sqrt{n}$ 점선 그래프를 추가하였습니다. 이항 분포의 성질을 이용해서 해석적인 방식으로 해를 구한 결과와 주방정식을 풀어서 확률 밀도 함수 이용하여 구한 결과는 당연히 일치합니다. 

이번에는 초기 확률 밀도 함수가 공간상의 두 곳에서 국소화 된 경우를 생각했습니다. 초기 조건은 $p(0, 0, 0) = \frac{}{2}, p(-4, -4, 0) = \frac{1}{2}$ 입니다. 원점과 $(-4, -4)$에 절반씩 확률 밀도 함수가 나누어진 상태에서 시간에 따른 변화 입니다. 주방정식의 선형성에 의해서 원점의 확률 밀도 함수의 전파와 $(-4, -4)$점의 확률 밀도 함수의 전파의 독립적인 합이라고 봐도 됩니다. 

 

위에서와 마찬가지로 $\langle r^2 \rangle(n)$을 구해 볼 수 있습니다. 이 경우에는 

(위 그래프의 $y$축은 $\sqrt{\langle r^2 \rangle}$ 입니다. )

 

와 같이 주어집니다. $n=0$에서 평균은 $(-2,2)$이지만, $(0, 0), (-4, -4)$에 절반의 확률이 분포되어 있기 때문에, $\langle r^2 \rangle(0)$ 값이 $0$이 아닌 값이 됩니다. 

 

위 결과는 $p(0, 0, 0) = 1$의 초기 조건에서의 확률 밀도 함수의 결과를 이용해서도 계산이 가능한데, $p(x, y, n)_0$을 이전 결과에서 구한 $p(0, 0, 0) = 1$의 초기 조건에서의 확률 밀도 함수라 하겠습니다. 이 문제에서 초기 조건은 $p(0, 0, 0) = \frac{1}{2}, p(-4, -4, 0) = \frac{1}{2}$이고 앞에서 언급한 것 처럼 임의의 $n$에서 확률 밀도 함수는 $p(x, y, n) = \frac{1}{2}p(x, y, n)_0 + \frac{1}{2}p(x+4, y+4, n)_0$이 됩니다. $p(x, y, n)$의 평균은 $n$에 상관 없이 $(x_m, y_m) = (-2, 2)$이 되고, $\langle r^2 \rangle (n)$은

$$\langle r^2 \rangle (n) = \sum_{x, y} ((x-x_m)^2 + (y-ym)^2)(\frac{1}{2}p(x, y, n)_0 + \frac{1}{2}p(x+4, y+4, n)_0 )$$

가 됩니다.  $(x-x_m)^2 = x^2 - 2x x_m + x_m^2$이므로

$$\langle (x-x_m)^2 \rangle  = \langle x^2 \rangle -2x_m \langle x \rangle + x_m^2$$

이 성립합니다. 앙상블 평균을 $p(x, y, n)_0$에 대해서 적용하면, $\langle x \rangle$은 $0$이기 때문에 $\langle (x-x_m)^2 \rangle = \langle x^2 \rangle + x_m^2$ 이 됩니다. $\langle r^2 \rangle = n$이므로, 우리가 구하려고 하는 $\langle r^2 \rangle(n)$값은

$$\langle r^2 \rangle (n) = n + \alpha$$

가 됩니다. 여기서 $\alpha$라고 쓴 값은 초기 분포 $p(x, y, 0)$에 의해 정해지는 상수 값 입니다. 위 분포에 대응되는 $\alpha$값은 $8$입니다. 왜 $8$이 되는지는 각자 계산해 보길 바랍니다. 

 

위 결과를 일반화 하면, 임의의 $n=0$에서의 확률 밀도 함수 $p(x, y, 0)$하에서 랜덤 워크를 한다면 임의의 $n$에서 $\langle r^2 \rangle(n)$은 $n + c$ 와 같이 주어지게 됩니다. 여기서 $c$는 앞에서 $\alpha$라고 쓴 것과 같은 상수의 역할을 하는 숫자 입니다. "상수" 라는 것을 강조하기 위해서 $c$라고 썼습니다. 즉, 초기 확률 밀도 함수의 개형과는 상관 없이 $n$이 충분히 커지면 $\langle r^2 \rangle(n) \approx n$이 성립합니다. 

 

치우친 랜덤 워크의 주방정식

 

위에서의 랜덤 워크는 "치우치지 않은" 랜덤 워크 였습니다. 즉, 이동할 수 있는 모든 방향으로 확률 밀도 함수가 동일한 비율로 그 다음 스텝으로 전파 되는 것 입니다. 1차원 정수 격자 공간에서 $x$에 할당 된 확률 밀도 함수 $p(x,n)$은 그 다음 스텝에서 $x=x-1$로 $\frac{1}{2} p(x, n)$이 전파 되고, $x = x+1$로 $\frac{1}{2} p(x, n)$이 전파 되는 것 입니다. 치우친 랜덤 워크에서는 이동 할 수 있는 방향 마다 그 비율이 달라 집니다. 이 때 주방정식은

$$p(x, n+1) = \alpha p(x-1, n) + (1-\alpha) p(x+1,n)$$

과 같이 나타낼 수 있습니다. 일반적으로 $\alpha \ne \frac{1}{2}$ 가 됩니다. 만일 $\alpha \gt \frac{1}{2}$가 된다면, (예를들어서 $\alpha = 0.8$을 생각할 수 있습니다) $x = x-1$의 확률 밀도 함수 중, 80%가 $x = x$로 옮겨지게 되고, 20%가 $x = x-1$로 옮겨지게 됩니다. 따라서 확률 밀도 함수의 평균은 오른쪽($x$값이 증가하는 방향)으로 이동하게 됩니다. 

 

2차원 정수 격자 랜덤 워크의 예시로 들면, 오른쪽-왼쪽으로 전파할 때, 그 비율을 (예를들면) $80 : 20$으로 위-아래로 전파할 때, 그 비율을 $60 : 40$ 이라고 한다면, 이 때의 확률 밀도 함수의 변화 $p(x, y, n)$은 아래 애니메이션과 같게 됩니다. 

위 애니메이션에서 볼 수 있든, 하늘색 원의 전체적인 움직임(평균)은 오른쪽 위를 향하고 있습니다. 이 처럼 치우쳐진 랜덤 워크를 주방정식의 "숫자 조절"을 통해서도 구현할 수 있습니다. 

 

아래는 위 애니메이션을 구현하기 위한 파이썬 프로그램 코드 입니다. 좀 더 공부를 하고 싶으신 분은 아래 코드를 이용하여 이것 저것 시도해 볼 수 있습니다. 

 

import numpy as np
import matplotlib.pyplot as plt
from matplotlib import animation

def update(num, list_list_x, list_list_y, list_list_value):
    Z = np.array([list_list_x[num], list_list_y[num]]).T
    point.set_offsets(Z)
    point.set_sizes(np.array(list_list_value[num]).T)

def return_index(x,y,L):
    return 2*L*y + x

def return_x_y(index, L):
    y, x = divmod(index, 2*L)
    return x-L, y-L

L = 15
N = 2*L * 2*L
A = np.zeros([N, N])

for x in range(2*L):
    for y in range(2*L):
        x_l, y_l = x-1, y
        x_r, y_r = x+1, y
        x_t, y_t = x, y+1
        x_b, y_b = x, y-1

        a = return_index(x, y, L)

        if 0<=x_l < 2*L and 0<=y_l < 2* L:
           b = return_index(x_l, y_l, L)
           A[a,b] = 0.40
           A[b,a] = 0.10

        if 0<=x_t < 2*L and 0<=y_t < 2* L:
           b = return_index(x_t, y_t, L)
           A[a,b] = 0.2
           A[b,a] = 0.30

p = np.zeros(N)
a = return_index(L, L, L)
p[a] =1

list_list_x, list_list_y, list_list_value =[], [], []

list_list_x.append([0])
list_list_y.append([0])
list_list_value.append([1.0*2000])

S = 14
for iteration in range(1, S):
    p = np.matmul(A, p)
    list_x, list_y, list_value = [], [], []
    for index, pp in enumerate(p):
        if not pp==0.0:
            x,y = return_x_y(index, L)
            list_x.append(x)
            list_y.append(y)
            list_value.append(pp*2000)

    list_list_x.append(list_x)
    list_list_y.append(list_y)
    list_list_value.append(list_value)

fig = plt.figure(figsize=(8, 8))
ax = plt.axes()
ax.grid()

ax.set_xlim(-L, L)
ax.set_ylim(-L, L)
ax.set_xlabel('x')
ax.set_ylabel('y')

point = ax.scatter([], [], s=50, color = 'dodgerblue')

ani = animation.FuncAnimation(fig, update, S, fargs=(list_list_x, list_list_y, list_list_value), interval=500, blit=False)
#ani.save('2dim_2.gif', writer=animation.PillowWriter(fps=2))
plt.show()