でたぁっ 感動と失敗の備忘録

データ解析を担当することになったサラリーマンの備忘録

【Rメモ】マルコフ連鎖とPageRank

 マーケティングではSTP戦略という手法がよく使われます。マーケティングの第一人者フィリップ・コトラー氏が提唱した代表的な手法。効果的に市場を開拓する為に 1.市場を顧客のニーズごとにセグメント化 2.参入すべきセグメントを選定しターゲットを明確にする 3.セグメントの中で顧客に対するベネフィット(利益)を検討し自らのポジションを確立する といった流れ。この戦略のS(セグメント)を議論する時、私はたたき台としてマルコフ連鎖をよく使っています。シンプルな手法ですが変化を可視化して共有すると議論がもりあがります。マーケティング担当者の受けも上々。今まではExcelで分析していたがRで完結させたくなりコードを書いたのでφ(..)メモメモ。

Rでマルコフ連鎖

 R+マルコフ連鎖については、teramonagiさんの資料 がよく分かります。(Rコードも記載されてます)。先ほどのマーケティングの例なら、現在のセグメント構成比が初期状態で、そのセグメント間の遷移確率が分かれば次の状態(セグメント構成比)を計算できるといったものです。計算は隣接行列×状態の構成比といった行列の演算、ExcelだとMMULT関数で一発。当然Rも行列の掛け算は %*% で行うことが出来るので一発。あとはその演算の繰り返し。自分で書くつもりでしたが、先ほどのteramonagiさんの資料にその関数が記載されていたのでそれを拝借。以下実際のデータとRコード。


初期構成比のデータ(セグメントの構成比)

V1 0.080
V2 0.060
V3 0.300
V4 0.100
V5 0.100
V6 0.360


遷移確率のデータ(各セグメント間で遷移する比率 ※V2->V3の確率は0.3)

V1 V2 V3 V4 V5 V6
V1 0.700 0.360 0.070 0.040 0.100 0.040
V2 0.100 0.240 0.090 0.080 0.050 0.010
V3 0.080 0.300 0.700 0.400 0.180 0.070
V4 0.020 0.050 0.050 0.250 0.170 0.060
V5 0.030 0.030 0.020 0.130 0.270 0.130
V6 0.070 0.020 0.070 0.100 0.230 0.690
1.000 1.000 1.000 1.000 1.000 1.000

setwd("d:/R/マルコフ連鎖") # 作業ディレクトリの変更

d1 <- as.matrix(read.table("遷移確率.txt", header=F)) # 遷移確率をmatrix形式で読込む
d2 <- as.matrix(read.table("初期構成比.txt", header=F)) # 初期構成比をmatrix形式で読込む

#----------------------------------------------------------
# 参考:http://www.slideshare.net/teramonagi/ss-5190440 p17
#
# 数日後の天気確率計算関数
# weather.initial_:
# transition.matrix_:お天気遷移行列
# size_:何日間の遷移を計算するのか
TranslateWeather <- function(weather.initial_, transition.matrix_, size_){
  result <- matrix(NA, nrow = size_ + 1, ncol = nrow(weather.initial_)) # 格納先
  result[1,] <- t(weather.initial_) # 初期値を結果に格納
  weather.day <- weather.initial_ # 初期値
  for(day in 1:size_){
    weather.day <- transition.matrix_ %*% weather.day # マルコフ連鎖(行列演算)
    result[day + 1,] <- t(weather.day) # 結果を転置して格納
  }
  return(result)
}
#----------------------------------------------------------
res <- TranslateWeather(d2, d1, 100) # 引数:初期値、遷移行列、回数
head(res)

          [,1]       [,2]      [,3]       [,4]       [,5]      [,6]
[1,] 0.0800000 0.06000000 0.3000000 0.10000000 0.10000000 0.3600000
[2,] 0.1270000 0.06600000 0.3176000 0.08320000 0.09700000 0.3092000
[3,] 0.1602880 0.07172200 0.3246640 0.07756200 0.08934400 0.2764200
[4,] 0.1838417 0.07589820 0.3280606 0.07418924 0.08359412 0.2544162
[5,] 0.2004804 0.07878419 0.3296510 0.07189805 0.07964253 0.2395439
[6,] 0.2121961 0.08075424 0.3302923 0.07031774 0.07696189 0.2294777

ggplot2で可視化

 先ほどのマルコフ連鎖の結果を ggplot2パッケージを使い可視化。(ggplot2で書きたくなったので)

library(ggplot2)
library(reshape)

res <- as.data.frame(TranslateWeather(d2, d1, 10)) # 引数:初期値、遷移行列、回数
res2 <- cbind(t=as.numeric(rownames(res)), res) # 行番号を列として設ける
res3 <- melt(res2, id="t") # ggplot2用のデータ形式に変換
head(res2)

  t        V1         V2        V3         V4         V5        V6
1 1 0.0800000 0.06000000 0.3000000 0.10000000 0.10000000 0.3600000
2 2 0.1270000 0.06600000 0.3176000 0.08320000 0.09700000 0.3092000
3 3 0.1602880 0.07172200 0.3246640 0.07756200 0.08934400 0.2764200
4 4 0.1838417 0.07589820 0.3280606 0.07418924 0.08359412 0.2544162
5 5 0.2004804 0.07878419 0.3296510 0.07189805 0.07964253 0.2395439
6 6 0.2121961 0.08075424 0.3302923 0.07031774 0.07696189 0.2294777

head(res3)

  t variable     value
1 1       V1 0.0800000
2 2       V1 0.1270000
3 3       V1 0.1602880
4 4       V1 0.1838417
5 5       V1 0.2004804
6 6       V1 0.2121961

gg2 <- ggplot(res3, aes(x=reorder(t, t), y=value, fill=variable)) # 積上げ棒グラフ reorderで横軸順序を指定
gg2 <- gg2 + geom_bar(colour=I("black"), alpha=0.8, linetype=1, size=0.7) # ちょっとキレイにする
plot(gg2)

f:id:deta:20130519040721p:plain
結構きれい。この遷移確率が続くといった仮定は入りますが、初期状態では小さいV1と大きなV6は逆転するといった結果を知ることが出来ます。このような積上げ棒グラフを描きたい場合、x軸の順序を並べ替える必要がある場合がありますのでご注意ください。実際にハマってしまいました。

マルコフ連鎖PageRank

 マルコフ連鎖は時間tを無限大にすると収束します。収束値は連立方程式で解けます。(これは会社の先輩から教えてもらいました ※平衡方程式というらしい) ご存じの方もいると思いますが、その収束値はPageRankと同じ値になります。私はたまたま知ったが。。。(1年くらい前にネットワーク図のノードの評価について駒場の学生さんとお話しする機会がありPageRankが話題に。職場に戻り算出していたマルコフ連鎖の収束値とPageRankの相関をとったら一致して気づく) RではigraphパッケージでPageRankを計算できるので一致するか確認。

library(igraph)
library(reshape)

d3 <- d1 # 遷移確率をコピー
rownames(d3) <- colnames(d3) # rownamesの設定
d3 <- as.data.frame(t(d3)) # from to 形式にする為の転置
d3 <- cbind(from=rownames(d3), d3) # 行名を列として設ける
d3 <- melt(d3, id="from") # from to weight 形式に整形
head(d3)

  from variable value
1   V1       V1  0.70
2   V2       V1  0.36
3   V3       V1  0.07
4   V4       V1  0.04
5   V5       V1  0.10
6   V6       V1  0.04

g <- graph.data.frame(d3[1:2], directed=T) # graph形式に変換
E(g)$weight <- d3[[3]] # weightは遷移確率
pr <- page.rank(g, directed=T, damping=1.0) # dampingのデフォルトは0.85
round(pr$vector, 5) # PageRankの結果

     V1      V2      V3      V4      V5      V6 
0.23924 0.08494 0.32928 0.06684 0.07135 0.20835

res <- TranslateWeather(d2, d1, 100) # 100回計算させる
tail(round(res, 5)) # マルコフ連鎖の結果

          [,1]    [,2]    [,3]    [,4]    [,5]    [,6]
[96,]  0.23924 0.08494 0.32928 0.06684 0.07135 0.20835
[97,]  0.23924 0.08494 0.32928 0.06684 0.07135 0.20835
[98,]  0.23924 0.08494 0.32928 0.06684 0.07135 0.20835
[99,]  0.23924 0.08494 0.32928 0.06684 0.07135 0.20835
[100,] 0.23924 0.08494 0.32928 0.06684 0.07135 0.20835
[101,] 0.23924 0.08494 0.32928 0.06684 0.07135 0.20835

PageRankの結果とマルコフ連鎖の収束値が一致していることが確認できます。(100回計算しなくても既に収束しているようですが) また、PageRankにはdampingパラメータがあります。デフォルトは0.85ですが1にしないとマルコフ連鎖の収束値にはなりません。

参考
マルコフ連鎖モンテカルロ法入門-1
reorder 関数を使って横軸の順序の並べ替えを行う方法