【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)
結構きれい。この遷移確率が続くといった仮定は入りますが、初期状態では小さい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にしないとマルコフ連鎖の収束値にはなりません。