Summary

본 내용은 STAT-110 강의 31 번의 일부 내용과 선형대수학(길버트 스트랭교수님) 교재의 일부 내용을 바탕으로 마코프 체인의 응용 예제를 작성합니다.

웹사이트 A, B, C 가 있다. 이 웹사이트의 10년 후 사용자수 예측과 50년 후 사용자수 예측을 마코프 체인과 전이행렬을 사용해서 예측을 수행한다.

이를 위해 아래 작업을 수행한다.

마코프 체인 모델 설정

웹사이트 A, B, C 가 있다. 매년 각 사이트의 사용자 이동수가 아래와 같다.

초기 사용자수는 \(u = [10000, 15000, 12000]^T\) 이다.

전이행렬 (Transition matrix)

전이 행렬을 Q 라고 할 경우, STAT-110 에서는 가로 행의 합이 1 이 되게 구성한다. 하지만, 여기서는 컬럼 기준으로 컬럼의 합이 1이 되게 구성한다.

\(Q = \begin{bmatrix} 0.7 & 0.15 & 0.05 \\ 0.2 & 0.7 & 0.1 \\ 0.1 & 0.15 & 0.85 \end{bmatrix}\)

set.seed(2021)

# 전이행렬 구성
Q <- matrix(c(0.7, 0.15, 0.05, 0.2, 0.7, 0.1, 0.1, 0.15, 0.85), nrow = 3, byrow = T)
Q # 전이행렬
##      [,1] [,2] [,3]
## [1,]  0.7 0.15 0.05
## [2,]  0.2 0.70 0.10
## [3,]  0.1 0.15 0.85
# 사용자수 구성. 만약 s 를 확률벡터로 구성하면, 단계가 진행할 수로 정상분포로 수렴함을 볼 수 있을 것이다.
s <- matrix(c(30000, 12000, 9000), ncol = 1) # 실제 사용자수로 구성함
#s <- matrix(c(0.5, 0.3, 0.2), ncol = 1) # 확률벡터로 구성함
s # 사용자 수
##       [,1]
## [1,] 30000
## [2,] 12000
## [3,]  9000

전이행렬 계산을 이용한 10년, 50년 후 사용자 수 구하기

# 전이행렬, 사용자수 벡터, 예측 기간
getUsers <- function(Q, s, num) {
  if (num == 1) {
    return (Q %*% s)
  }
  
  for(i in 1:num) {
    s <- Q %*% s
  }
  return (s)
}

# 1 년 뒤 사용자수
#s1 <- Q %*% s
s_result <- getUsers(Q, s, 1)
s_result # 1년뒤 사용자 수
##       [,1]
## [1,] 23250
## [2,] 15300
## [3,] 12450
# 10 년 뒤 사용자수
s_result <- getUsers(Q, s, 10)
s_result # 10년뒤 사용자 수
##          [,1]
## [1,] 12181.39
## [2,] 15987.23
## [3,] 22831.38
# 50 년 뒤 사용자수
start_time <- Sys.time()
s_result <- getUsers(Q, s, 50)
Sys.time() - start_time
## Time difference of 0.002243996 secs
s_result # 50년뒤 사용자 수
##          [,1]
## [1,] 11769.23
## [2,] 15692.31
## [3,] 23538.46

전이행렬의 고유값, 고유벡터

전이행렬 Q 는 몇 가지 특성을 가진다.

eig <- eigen(Q)
my_lambda <- eig$values
S <- eig$vectors

# 고유값 확인
my_lambda # 전이행렬의 고유값
## [1] 1.0000000 0.7280776 0.5219224
# 고육벡터 확인
S # 전이행렬의 고유벡터 
##           [,1]       [,2]       [,3]
## [1,] 0.3841106 -0.4573522  0.6064256
## [2,] 0.5121475 -0.3570898 -0.7766956
## [3,] 0.7682213  0.8144420  0.1702700

전이행렬의 대각화를 이용한 10년, 50년 후 사용자 계산

초기 사용자 수는 고유벡터 행렬의 선형 결합으로 표현이 가능하다.

\(\vec{u} = S\vec{c}\)

따라서, 계수 벡터는 \(\vec{c} = S^{-1}\vec{u}\) 로 구할 수 있다.

그리고 k 년 후의 값은 위에서 구한 전이행렬의 고유값, 고유벡터를 이용해서 아래와 같이 구한다.

\(u_k = A^k u_0 = S \Lambda^k S^{-1} u_0 = S \Lambda^k \vec{c} = c_1(\lambda_1)^kx_1 + c_2(\lambda_2)x_2 + c_3(\lambda_3)x_3\)

(Gilbert Strang, 2009, Introduction to Linear Algebra, Wellesley-Cambridge Press, 431-434)

# 계수벡터를 구한다.
c_vec <- solve(S) %*% s
c_vec # 초기값 s 를 이용해서 계수벡터를 구한다.
##           [,1]
## [1,]  30640.21
## [2,] -20848.61
## [3,]  14339.12
# 예측 년도 
getUsers_eig <- function(k) {
  u_k <- c_vec[1] * my_lambda[1]^k * S[,1] + 
    c_vec[2] * my_lambda[2]^k * S[,2] + 
    c_vec[3] * my_lambda[3]^k * S[,3]
}

# 10년 후
u_k <- getUsers_eig(10)
u_k # 10년후 사용자 수
## [1] 12181.39 15987.23 22831.38
# 50년 후
start_time <- Sys.time()
u_k <- getUsers_eig(50)
Sys.time() - start_time
## Time difference of 0.008543968 secs
u_k # 50년 후 사용자 수 
## [1] 11769.23 15692.31 23538.46
# 비율로 나타내면, 고유벡터 x1 과 같음을 확인한다.
s_result / sum(s_result) # 행렬 계산을 통해서 구한 사용자수 확률벡터
##           [,1]
## [1,] 0.2307693
## [2,] 0.3076923
## [3,] 0.4615384
u_k / sum(u_k) # 고유값, 고유벡터 계산을 통해서 구한 사용자수 확률벡터
## [1] 0.2307693 0.3076923 0.4615384
S[,1] / sum(S[,1]) # 전이행렬 Q에서 람다=1을 가진 고유벡터의 확률벡터
## [1] 0.2307692 0.3076923 0.4615385
# 위에서 구한 확률벡터가 정상분포임을 확인한다. Q*s = s
sp <- matrix(c(0.2307692, 0.3076923, 0.4615385), ncol = 1)
Q %*% sp # 정상분포를 전이행렬과 곱하면 정상분포가 나온다.
##           [,1]
## [1,] 0.2307692
## [2,] 0.3076923
## [3,] 0.4615385

오래 실행할 경우, Steady state 에 접근한다.

행렬을 직접 곱하는 방식으로, 오래 실행한다. 그리고 이 결과로부터 각 사이트를 사용하는 사용자 수의 분포가 일정하게 유지됨을 확인한다. 또한 사이트 C 가 사용자수가 많으므로 랭킹이 높다는 것을 알 수 있다.

# 50 년간 각 사이트의 사용자수 추이
A_users <- c()
B_users <- c()
C_users <- c()
years <- 50
for(i in 1:years) {
  s_result <- getUsers(Q, s, i)
  A_users <- c(A_users, s_result[1])
  B_users <- c(B_users, s_result[2])
  C_users <- c(C_users, s_result[3])
}

# 추이를 보기 위한 데이터 셋 구축
df_users <- data.frame(years = c(rep(1:years, 3)), 
      users = c(A_users, B_users, C_users), 
      site = c(rep("A", years), rep("B", years), rep("C", years)))

#df_users

library(ggplot2)
g <- ggplot(df_users, aes(x = years, y = users, colour = site)) + geom_line()
g

simulation 을 통한 사용자 분포

마코프 체인에 따라서, 사용자의 현재 사이트가 주어졌을 때, 다음 사이트를 구한다. 그리고 방문한 사이트의 hit 수를 저장한다. 이렇게 오랫동안 실행했을 경우, 각 사이트의 방문수를 막대 그래프로 확인한다.

그리고 이 결과가 위에서 구한 값과 비슷함을 확인한다.

# 사이트를 입력받으면, 다음 이동할 사이트를 반환한다.
next_site <- function(current_site) {
  if(current_site == 'A') {
    # 만약 현재 상태가 A 라면, (A, B, C) 중 하나를 선택하는데
    # A 일때의 전이확률 Q[,1] 에 따른다.
    n_site <- sample(c('A', 'B', 'C'), 1, prob = Q[,1])
  } else if (current_site == 'B') {
    n_site <- sample(c('A', 'B', 'C'), 1, prob = Q[,2])
  } else {
    n_site <- sample(c('A', 'B', 'C'), 1, prob = Q[,3])
  }
}

# 각 사이트 방문 횟수를 저장할 named vector 를 0으로 초기화한다.
user_count <- c('A' = 0, 'B' = 0, 'C' = 0)
current_site <- 'A' # A 부터 시작한다.
for (i in 1:100000) {
  user_count[current_site] <- user_count[current_site] + 1
  current_site <- next_site(current_site)
}
#user_count / sum(user_count)
barplot(user_count)

REFERENCES

[1] Gilbert Strang, 2009, Introduction to Linear Algebra, Wellesley-Cambridge Press, 431-434