某商圈有 3 家店均販賣 3C 產品，但各自規模及經營方式有所差異，市場調查專家欲從客戶交易行為了解各店家特色是否有效，記錄某月客戶交易行為如下：

店一： 前月客戶數=350 月底客戶數=380 新增客戶中： 來自店二有 40 來自店三有 30 流失客戶中： 去店二有 30 去店三有 10

店二： 前月客戶數=400 月底客戶數=390 新增客戶中： 來自店一有 30 來自店三有 20 流失客戶中： 去店一有 40 去店三有 20

店三： 前月客戶數=250 月底客戶數=230 新增客戶中： 來自店一有 10 來自店二有 20 流失客戶中： 去店一有 30 去店二有 20

(1) 假使各店經營型態不變，則此商圈在未來 10 個月各店客戶佔有率如何？

(2) 此商圈各店知客戶群是否有固定行為模式？若有，最終分別在各店客戶佔有率為何？

總客戶數 350 + 400 + 250 = 1000 店一 店二 店三 前月客戶佔有率 350/1000 400/1000 250/1000 忠誠客戶數 350 - (30 + 10) = 310 400 - (40 + 20) = 340 250 - (30 + 20) = 200 前月客戶數比率 310/350 340/400 200/250

客戶移轉率 店一 店二 店三 店一 310/350 30/350 10/350 店二 40/400 340/400 20/400 店三 30/250 20/250 200/250

---

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
}

NumOfCust <- c(350, 400, 250)
TotalCust <- sum(NumOfCust)
CustMoveA2BC <- c(30, 10)
CustMoveB2AC <- c(40, 20)
CustMoveC2AB <- c(30, 20)

P <- matrix(c(
      (NumOfCust[1]-sum(CustMoveA2BC))/NumOfCust[1], 
              CustMoveA2BC[1]/NumOfCust[1], 
              CustMoveA2BC[2]/NumOfCust[1], 
              CustMoveB2AC[1]/NumOfCust[2], 
      (NumOfCust[2]-sum(CustMoveB2AC))/NumOfCust[2], 
              CustMoveB2AC[2]/NumOfCust[2], 
              CustMoveC2AB[1]/NumOfCust[3], 
              CustMoveC2AB[2]/NumOfCust[3], 
      (NumOfCust[3]-sum(CustMoveC2AB))/NumOfCust[3]), 
            byrow=TRUE, nrow=3, ncol=3)

P0 <- matrix(c(NumOfCust[1]/TotalCust, NumOfCust[2]/TotalCust, NumOfCust[3]/TotalCust), 
             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", i, ans[1], ans[2], ans[3]))

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

01) 店一:0.3800 店二:0.3900 店三:0.2300
02) 店一:0.4032 店二:0.3825 店三:0.2144
03) 店一:0.4211 店二:0.3768 店三:0.2021
04) 店一:0.4349 店二:0.3725 店三:0.1926
05) 店一:0.4455 店二:0.3693 店三:0.1851
06) 店一:0.4538 店二:0.3669 店三:0.1793
07) 店一:0.4601 店二:0.3651 店三:0.1747
08) 店一:0.4650 店二:0.3638 店三:0.1712
09) 店一:0.4688 店二:0.3628 店三:0.1684
10) 店一:0.4717 店二:0.3620 店三:0.1663

<穩態> 店一:0.4815 店二:0.3598 店三:0.1587
<穩態需要移轉次數> 店一:2.0768 店二:2.7793 店三:6.3012

由右圖脊線走向可知，移轉機率朝向穩態發展，但需要超過 10 個月以上才會達到穩態。