科學家探索失智老鼠行為模式，設計 3 間不同密室各自擺上不同氣味食物，觀察這隻失智老鼠在各密室活動情形，並加以記錄分析發現以下模式：

(1) 此失智老鼠自【密室一】出發，在 10 步驟內分別停留在各密室機率如何？

(2) 此失智老鼠是否有固定行為模式？若有，最終分別在各密室機率為何？

options(digits=4)

rm(list=ls())

Pn <- function(P, n) {
  ans <- P
  
  if (n > 1) {
    for (i in seq(2, n)) {
      # 矩陣相乘
      ans <- P %*% ans
    }
  }
  
  ans
}

P <- matrix(c(1/2, 1/4, 1/4, 0, 1/3, 2/3, 1/2, 1/2, 0), byrow=TRUE, nrow=3, ncol=3)
P0 <- matrix(c(1, 0, 0), byrow=TRUE, nrow=1, ncol=3)

# 移轉機率結果矩陣
A <- matrix(0, nrow=10, ncol=3)

for (i in 1:10) {
  # 矩陣相乘
  ans <- P0 %*% Pn(P, i)
  A[i,] <- c(ans[1], ans[2], ans[3])
  print(sprintf("%02d) 密室一:%.4f 密室二:%.4f 密室三:%.4f", i, ans[1], ans[2], ans[3]))
}

# 移轉機率結果矩陣欄位名稱
colnames(A) <- c('密室一', '密室二', '密室三')

# 繪製移轉機率結果矩陣立體觀點圖
persp(A, xlab='密室一', ylab='密室二', zlab='密室三', main='機率趨勢平面圖', 
      theta=30, phi=30, expand=0.5, col="lightblue",
      ltheta= 120, shade=0.75, ticktype="detailed")

# 狀態關係： 轉置矩陣 - 單位矩陣
W <- t(P) - diag(nrow=3, ncol=3)

# 各狀態機率和 = 1
W2 <- rbind(W, c(1, 1, 1))

# 解穩態方程式： 移除首列多餘關係式
ans <- solve(W2[2:4, ], c(0, 0, 1))

#穩態機率
print(sprintf("<穩態> 密室一:%.4f 密室二:%.4f 密室三:%.4f", ans[1], ans[2], ans[3]))

#期望再現時間
print(sprintf("<穩態需要移轉次數> 密室一:%.4f 密室二:%.4f 密室三:%.4f", 1/ans[1], 1/ans[2], 1/ans[3]))

01) 密室一:0.5000 密室二:0.2500 密室三:0.2500
02) 密室一:0.3750 密室二:0.3333 密室三:0.2917
03) 密室一:0.3333 密室二:0.3507 密室三:0.3160
04) 密室一:0.3247 密室二:0.3582 密室三:0.3171
05) 密室一:0.3209 密室二:0.3591 密室三:0.3200
06) 密室一:0.3204 密室二:0.3599 密室三:0.3196
07) 密室一:0.3200 密室二:0.3599 密室三:0.3201
08) 密室一:0.3200 密室二:0.3600 密室三:0.3199
09) 密室一:0.3200 密室二:0.3600 密室三:0.3200
10) 密室一:0.3200 密室二:0.3600 密室三:0.3200

<穩態> 密室一:0.3200 密室二:0.3600 密室三:0.3200
<穩態需要移轉次數> 密室一:3.125 密室二:2.7778 密室三:3.125

由右圖脊線走向可知，移轉機率朝向穩態發展。