07wk-1: 마코프체인 (2)

최규빈

2023-04-13

강의영상

youtube: https://youtube.com/playlist?list=PLQqh36zP38-wi47mCKi03xoqvwIzkxG0H

import

import numpy as np

확률과정

- 동전을 무한히 던지는 시행을 생각하자. 동전을 10번 던져서 결과를 관찰했다고 하자. 동전을 30번째 던져서 앞면이 나올지 뒷면이 나올지 알고 싶다면?

- 현재 삼성전자 주가는 66000이다. 20일뒤의 삼성전자 주가가 얼마일지 알고 싶다면?

- 원래 미래를 예측하기 위해서 해야하는 과정

그림1: 1400만개의 미래를 탐색중인 Doctor Strange

- 하지만 현실적으로는 이게 너무 힘들지 않을까?

날씨예측

- 아래와 같이 세상의 법칙이 있다고 하자.

• 어제 맑음 $\to$ 오늘도 맑음: 40% // 오늘은 비: 60%
• 어제 비 $\to$ 오늘은 맑음: 70% // 오늘도 비 30%

- 모든 $t$에 대하여 확률변수 $X_t$를 아래와 같이 정의하자.

• $X_t=\{\begin{array}{ll}0& \text{맑음}\\ 1& \text{비}\end{array}$

- 오늘 (2023년4월13일) 비가 왔다고 치자. 10000일 뒤에도 비가 올 확률은 얼마일까?

풀이1

- $X_t=0$ 이라면? ($t$시점에 비가 오지 않았다면?)

np.random.rand() < 0.6 

False

- $X_t=1$ 이라면? ($t$시점에 비가 왔다면?)

np.random.rand(0) < 0.3

array([], dtype=bool)

- 두 코드를 합쳐보자.

def rain(before):
    if before == True: # 비가 왔음 
        after = np.random.rand() < 0.3
    else: # 비가 안왔음 
        after = np.random.rand() < 0.6 
    return after 

- 테스트

# 비가 왔음, Xt = 1 
sum([rain(1) for i in range(100)])

30

# 비가 안왔음, Xt = 0 
sum([rain(0) for i in range(100)])

60

- 하나의 $\omega$에 대응하는 길이가 10000인 확률과정을 관찰

def doctor_strange(today):
    lst = [today]
    for i in range(10000): 
        lst.append(rain(lst[i]))
    return lst 

today = True # 오늘 비가 왔다는 뜻 
arr = doctor_strange(today)

len(arr)

10001

- 4305개의 $\omega$에 대응하는 길이가 10000인 확률과정을 관찰

today = True # 오늘 비가 왔다는 뜻 
arr = np.array([doctor_strange(today) for ω in range(4305)])

arr[:,-1].mean()

0.4662020905923345

- 10000일 뒤에도 비가 올 확률은 약 46% 정도 인듯

풀이2

- 세상의 법칙을 다시 정리해보자.

• $X_{t-1}=0\Rightarrow X_t\sim Ber(0.6)$
• $X_{t-1}=1\Rightarrow X_t\sim Ber(0.3)$

- 정리하면

• $P(X_t=0)=P(X_{t-1}=0)\times 0.4+P(X_{t-1}=1)\times 0.7$
• $P(X_t=1)=P(X_{t-1}=0)\times 0.6+P(X_{t-1}=1)\times 0.3$

- 매트릭스형태로 표현하면

• $\left[\begin{array}{c}P(X_t=0)\\ P(X_t=1)\end{array}\right]=\left[\begin{array}{cc}0.4& 0.7\\ 0.6& 0.3\end{array}\right]\left[\begin{array}{c}P(X_{t-1}=0)\\ P(X_{t-1}=1)\end{array}\right]$
• ${\mathit{\mu }}_t=\mathbf{P}{\mathit{\mu }}_{t-1}$

- 이렇게 놓고 보니까

• ${\mathit{\mu }}_1=\mathbf{P}{\mathit{\mu }}_0$
• ${\mathit{\mu }}_2=\mathbf{P}{\mathit{\mu }}_1={\mathbf{P}}^2{\mathit{\mu }}_0$
• $\dots$
• ${\mathit{\mu }}_{10000}={\mathbf{P}}^{10000}{\mathit{\mu }}_0$

- 이제 계산을 해보자.

μ0 = np.array([[0],[1]])
μ0

array([[0],
       [1]])

P = np.array([[0.4,0.7],[0.6,0.3]])
P

array([[0.4, 0.7],
       [0.6, 0.3]])

P@P # P의 제곱

array([[0.58, 0.49],
       [0.42, 0.51]])

P@P@P@P # P의 4제곱

array([[0.5422, 0.5341],
       [0.4578, 0.4659]])

P@P@P@P @ P@P@P@P # P의 8제곱 

array([[0.53849182, 0.53842621],
       [0.46150818, 0.46157379]])

P@P@P@P@P@P@P@P @ P@P@P@P@P@P@P@P # P의 16제곱 

array([[0.53846154, 0.53846154],
       [0.46153846, 0.46153846]])

$\mathbf{P}$가 수렴하는거 같지 않어?

P@P@P@P@P@P@P@P@P@P@P@P@P@P@P@P @ P@P@P@P@P@P@P@P@P@P@P@P@P@P@P@P 

array([[0.53846154, 0.53846154],
       [0.46153846, 0.46153846]])

대충 ${\mathbf{P}}^{10000}\approx {\mathbf{P}}^{32}$

Plim = P@P@P@P@P@P@P@P@P@P@P@P@P@P@P@P @ P@P@P@P@P@P@P@P@P@P@P@P@P@P@P@P 
Plim 

array([[0.53846154, 0.53846154],
       [0.46153846, 0.46153846]])

Plim @ μ0

array([[0.53846154],
       [0.46153846]])

- 이 풀이에 따르면 10000일 뒤에 비가 올 확률은 46% 정도이다.

풀이3

- 세상의 법칙을 다시 정리해보자.

• $X_{t-1}=0\Rightarrow X_t\sim Ber(0.6)$
• $X_{t-1}=1\Rightarrow X_t\sim Ber(0.3)$

- 추측: 10000일 뒤에 비가 올 확률이 $p$라고 치자. 그렇다면 9999일 뒤에 비가 올 확률도 $p$ 아닐까?

이걸 가정하고 계산해보자

1. 9999일 뒤에 비가 안 올 확률 $1-p$

• 9999일 뒤에 비가 안오고, 10000일 뒤에는 비가 올 확률: $0.6(1-p)$
• 9999일 뒤에 비가 안오고, 10000일 뒤에는 비가 안 올 확률: $0.4(1-p)$

2. 9999일 뒤에 비가 올 확률 $p$

• 9999일 뒤에 비가 오고, 10000일 뒤에도 비가 올 확률: $0.3p$
• 9999일 뒤에 비가 오고, 10000일 뒤에는 비가 안 올 확률: $0.7p$

따라서 $0.6(1-p)+0.3p=p$

풀어보면 $0.6/1.3=p$

0.6/1.3

0.4615384615384615

풀이4

np.mean(doctor_strange(True)[1:])

0.462