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{cases} 0 & \text{맑음} \\ 1 & \text{비} \end{cases}\)

- 오늘 (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\)

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

  • \(\begin{bmatrix} P(X_t=0) \\ P(X_t=1) \end{bmatrix}= \begin{bmatrix} 0.4 & 0.7 \\ 0.6 & 0.3 \end{bmatrix} \begin{bmatrix} P(X_{t-1}=0) \\ P(X_{t-1}=1) \end{bmatrix}\)
  • \({\boldsymbol \mu}_t = {\bf P} {\boldsymbol \mu}_{t-1}\)

- 이렇게 놓고 보니까

  • \({\boldsymbol \mu}_1 ={\bf P}{\boldsymbol \mu}_0\)
  • \({\boldsymbol \mu}_2 ={\bf P}{\boldsymbol \mu}_1={\bf P}^2{\boldsymbol \mu}_0\)
  • \(\dots\)
  • \({\boldsymbol \mu}_{10000} ={\bf P}^{10000}{\boldsymbol \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]])

\({\bf 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]])

대충 \({\bf P}^{10000} \approx {\bf 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