내용
마코브 체인(Markov Chains)
Andrey Markov의 이름을 딴 Markov 체인은 한 "상태"(상황 또는 값 집합)에서 다른 "상태"로 이동하는 수학적 시스템입니다. 예를 들어, 아기의 행동에 대한 Markov 체인 모델을 만든 경우 "놀기", "먹기", "자고 있음" 및 "울음"을 상태로 포함할 수 있으며 다른 행동과 함께 상태 공간(state space)을 형성할 수 있습니다. 또한, 상태 공간의 맨 위에 있는 Markov 체인은 한 상태에서 다른 상태로의 "전환" 확률을 알려줍니다. 예를 들어 현재 놀이 중인 아기가 다음 상태에서 먼저 울지 않고 5분이내에 잠들 확률을 계산할 수 있습니다.
위 그림은 두 개의 상태(A 및 B)가 있는 경우로 4개의 가능한 전환이 있습니다(상태가 다시 자체로 전환될 수 있기 때문에 2가 아님). 우리가 'A'에 있으면 'B'로 전환하거나 'A'에 머무를 수 있습니다. 'B'에 있으면 'A'로 전환하거나 'B'에 머무를 수 있습니다. 이 두 상태 다이어그램에서 어떤 상태에서 다른 상태로 전환될 확률은 0.5입니다.
상태의 전환을 나타내기 위해 위의 Markov 체인 다이어그램 대신에 다음과 같은 전환행렬을 사용합니다. 이 행렬은 상태 전환의 확률을 즉, 전환확률을 계산하기 위해 사용합니다.
| A | B | |
|---|---|---|
| A | P(A|A):0.5 | P(B|A):0.5 |
| B | P(A|B):0.5 | P(B|B):0.5 |
이와 같은 마코브 체인은 생물학, 비즈니스, 화학, 공학, 물리학 등 다양한 상황에 적용되는 수학적 모델로서 동일한 방식으로 여러 번 수행되는 실험 또는 측정을 설명하는 데 사용됩니다. 여기서 각 시도(trial)의 결과는 가능한 여러 결과 중 하나가 되고 한 번의 시도의 결과는 단지 직전 시도에만 의존된다는 조건 하에 사용됩니다. 몇 가지 예를 통해 마코프 모델에 대해 알아봅니다.
예)
한 지역의 흡연 비율은 다음과 같이 나타낼 수 있습니다.
import numpy as np import numpy.linalg as la import pandas as pd import sympy as sp
x0=np.array([[0.61],[0.39]]);x0
array([[0.61],
[0.39]])
총 인구 100%에서 각 항목의 비율을 나타낸 것으로 위 열벡터의 합은 1이 됩니다. 이렇게 합이 1이 되는 벡터를 확률 벡터(probability vector)라고 합니다. 이러한 확률 벡터들로 구성된 행렬은 확률 행렬(stochastic matrix) 또는 전이 행렬(transition matrix)라고 합니다.
- 각 값은 [0, 1] 범위에서 존재
- 각 열벡터의 합은 1
마코브 체인은 확률행렬 P와 함께 고려되는 확률벡터 x0, x1, …등의 연속적인 수열 입니다.
$$\begin{align} x_1 &= Px_0\\ x_2 &= Px_1\\ x_3 &= Px_2\\ &\vdots\\ → &x_{k+1} = Px_k k=0, 1, 2, \cdots \end{align}$$$\mathbb{R}^n$ 벡터들의 마코브 체인이 어떠한 시스템 또는 일련의 실험을 설명할 때 xk의 각 항목은 그 시스템이 n개의 가능한 상태 내에 존재할 확률 또는 그 실험의 결과가 n개의 가능한 결과들 중의 하나일 확률을 나열합니다. 이러한 이유로 xk는 상태 벡터(state vector)라고 합니다.
매년 흡연과 금연의 평균 인구 비율의 변화는 다음과 같이 나타낼 수 있습니다.
| To | From | 흡연 | 금연 |
|---|---|---|
| 흡연 | 0.97 | 0.05 |
| 금연 | 0.03 | 0.95 |
| 합 | 1 | 1 |
위 데이터는 각 열의 합이 1이 되므로 다음과 같이 확률행렬 P로 고려할 수 있습니다.
P=np.array([[0.97, 0.05],[0.03, 0.95]]); P
array([[0.97, 0.05],
[0.03, 0.95]])
현재 흡연과 금연의 비율이 위 확률 벡터 x0와 같다면 다음 해의 비율은 마코브 체인에 의해 다음과 같이 계산됩니다.
x1=np.dot(P, x0); x1
array([[0.6112],
[0.3888]])
위 결과는 흡연은 61.12%, 금연은 38.88%가 될 것으로 예측합니다. 위와 동일한 계산에 의해 그 다음 년도의 비율을 계산할 수 있습니다.
위 계산을 반복하면 어떤 값에 수렴함을 알 수 있습니다.
re, P
(array([[0.61],
[0.39]]),
array([[0.97, 0.05],
[0.03, 0.95]]))
re={0:x0}
for i in range(1, 201):
re[i]=np.dot(P, re[i-1])
if i % 20 ==0:
print(f'{i}: {re[i]}')
20: [[0.6221696]
[0.3778304]]
40: [[0.62446592]
[0.37553408]]
60: [[0.62489922]
[0.37510078]]
80: [[0.62498098]
[0.37501902]]
100: [[0.62499641]
[0.37500359]]
120: [[0.62499932]
[0.37500068]]
140: [[0.62499987]
[0.37500013]]
160: [[0.62499998]
[0.37500002]]
180: [[0.625]
[0.375]]
200: [[0.625]
[0.375]]
위 결과는 확률행렬과 이전확률과의 행렬곱이 이전확률과 같게되는 부분으로 각 확률(흡연 62.5%, 금연 37.5%)이 수렴됨을 나타냅니다. 이와 같이 수렴되는 확률 벡터를 정상상태 벡터(steady-sate vector 또는 equilibrium vector)라고 합니다. 모든 확률 행렬은 정상 상태 벡터를 가질 수 있으며 식 1로 표현됩니다.
$$\begin{align}\tag{1} &Px = x\\ &Px-Ix = 0\\ &(P-I)x = 0 \end{align}$$식 1은 고유값과 고유벡터로 생성할 수 있는 특성방정식입니다. 즉, P의 고유값과 그에 대응하는 고유벡터가 이 방정식의 해로서 정상벡터를 계산할 수 있습니다. 확률행렬 P의 고유값 1에 대응하는 고유벡터를 찾을 수 있습니다. 위 확률행렬 P의 고유값과 고유행렬은 다음과 같습니다.
eigVal, eigVec=la.eig(P) eigVal
array([1. , 0.92])
np.around(eigVec, 3)
array([[ 0.857, -0.707],
[ 0.514, 0.707]])
고유값 1에 대응하는 고유벡터는 eigVec[:,1]이며 확률벡터로 수정하면 다음과 같이 위의 결과와 같습니다.
(eigVec[:,0]/np.sum(eigVec[:,0])).reshape(-1,1)
array([[0.625],
[0.375]])
예)
다음 중 확률행렬을 결정하여 봅니다.
확률행렬의 조건에 부합하는 행렬은 B입니다.
예)
다음 행렬 p의 정상상태 벡터를 계산합니다.
p=np.array([[0.2, 0.6], [0.8, 0.4]]) eigVal, eigVec=la.eig(p) eigVal
array([-0.4, 1. ])
np.around(eigVec, 3)
array([[-0.707, -0.6 ],
[ 0.707, -0.8 ]])
q=eigVec[:,1]/np.sum(eigVec[:,1]) np.around(q.reshape(2,1), 3)
array([[0.429],
[0.571]])
예)
의회선거에서 D, R, L당에 투표하는 비율에 대한 것으로 결과는 직전 결과에만 의존한다고 가정하면 이 일련의 결과는 마코브체인이 될 수 있습니다. 일반적으로 투표 경향은 다음 표와 같다고 가정하면 이 체인의 확률 행렬의 다음과 같습니다.
| To/From | D | R | L |
|---|---|---|---|
| D | 0.8 | 0.2 | 0.3 |
| R | 0.1 | 0.7 | 0.3 |
| L | 0.1 | 0.1 | 0.4 |
P=np.array([[0.8, 0.2, 0.3],[0.1, 0.7, 0.3],[0.1, 0.1, 0.4]]);P
array([[0.8, 0.2, 0.3],
[0.1, 0.7, 0.3],
[0.1, 0.1, 0.4]])
P와 같은 전이 행렬을 상수로 고려할 수 있다면 투표 결과인 일련의 벡터들은 마코브 체인을 따릅니다. 한 선거의 결과가 x0과 같다면, 다음(x1)과 그 다음(x2)의 결과?
x0=np.array([0.55, 0.40, 0.05]).reshape(3,1); x0
array([[0.55],
[0.4 ],
[0.05]])
x1=np.dot(P, x0);x1
array([[0.535],
[0.35 ],
[0.115]])
x2=np.dot(P, x1);x2
array([[0.5325],
[0.333 ],
[0.1345]])
위 시스템으로부터 정상상태 벡터(q)를 결정하여 봅니다.
eigVal, eigVec=la.eig(P) eigVal
array([1. , 0.6, 0.3])
np.around(eigVec, 3)
array([[-0.836, -0.707, -0.267],
[-0.502, 0.707, -0.535],
[-0.223, -0. , 0.802]])
q=eigVec[:,0]/np.sum(eigVec[:,0]) np.around(q.reshape(3,1), 3)
array([[0.536],
[0.321],
[0.143]])
np.sum(q)
1.0
정방행렬 T의 거듭제곱(Tn)의 모든 원소가 양수일 경우 그 행렬을 regular라고 합니다. regular인 확률행렬에 의해 나타낼 수 있는 마코브 체인을 regular Markov chain이라고 합니다.
예)
다음 행렬들 중에 regular 확률행렬을 결정해봅니다.
다음은 행렬 A에 대한 것으로 그 행렬 자체가 확률 행렬이고 모든 요소가 양수이므로 regular입니다. 행렬 B의 경우는 확률 행렬이지만 0을 포함하므로 그 자체가 regular는 아닙니다. 그러므로 거듭제곱의 형태를 확인해보면 다음과 같이 제곱일 경우 모든 원소들이 양수인 확률 행렬이 됩니다. 즉, regular입니다.
행렬의 거듭제곱은 numpy.linalg 모듈의 matrix_power()함수로 계산할 수 있으며 객체내의 원소를 비교하여 True/False를 반환하는 numpy 함수 any()를 사용하여 regular를 판단할 수 있는 함수를 생성하여 적용합니다.
B=np.array([[0.0, 0.2],[1.0, 0.8]]) la.matrix_power(B, 2)
array([[0.2 , 0.16],
[0.8 , 0.84]])
la.matrix_power(B, 3)
array([[0.16 , 0.168],
[0.84 , 0.832]])
def Regular(M, n):
re=[]
for i in range(2, n):
x=la.matrix_power(M, i)
if np.any(x<0) or np.any(x == 0):
re.append(i)
if len(re)==0:
print('regular')
else:
print("not regular")
print(re)
Regular(B, 101)
regular
위와 같이 행렬 C와 D는 제곱부터 0을 포함하므로 regular가 아닙니다.
C=np.array([[0.8, 0.2],[0, 1]])
D=np.array([[1, 0.1, 0.3],[0, 0.2, 0.5],[0, 0.7, 0.2]])
Regular(C, 101)
not regular
[2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99, 100]
Regular(D, 101)
not regular
[2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99, 100]
la.matrix_power(D, 20)
array([[1. , 0.98980775, 0.99138598],
[0. , 0.00466846, 0.00394556],
[0. , 0.00552379, 0.00466846]])
댓글
댓글 쓰기