2/6 に 「パターン認識と機械学習」(PRML)読書会 #11 ＠サイボウズ・ラボに のこのこ行ってきました。お疲れ様＞各位
今回は８章「グラフィカルモデル」の後半＋９章の K-means まで。
sum-product(積和アルゴリズム) や max-sum で、グラフィカルモデルが周辺化や同時分布の大域最大解を求めるためのツールになる、というお話。
いつものように、気付いたことその他箇条書きで。
- 有向グラフを因子グラフに変換するために必要なモラル化の説明で、sleepy_yoshi さんの資料に超ナツカシ漫画が出てきて、ずっこけたｗ 断じてアメフト漫画ではありませんｗｗｗ ちなみに、件の漫画は男子校の友人から貸してもらって読んだ……いやいや、もちろん純粋に内容が面白いからであって、そんな勘ぐ(ry
- sum-product でも max-sum でも、その導出のキーとなっているのが「積和交換」(例えば 8.63 式)なのだが、その操作(ΣとΠの順序を入れ替え。もちろん一般にそんなことは成立しない)については 8.52〜8.53 式で流すような感じで説明しているだけなので、読む人がおおむね混乱に陥っているようだった。両アルゴリズムのキーなのでもう少し丁寧に書いてあっても良かったんじゃあないかなあ。
- 8.65 式で「 G_m(x_m, X_m) ってなんだ！？」という議論に。因子グラフは「こういう因子分解がある」ことしか言っておらず、その因子がどういう形なのかは全く規定していない。つまり因子群 G_m は存在はわかるが、具体的な内容は不明*1。ただ一つだけ間違いないことは、G_m が先の「積和交換」が可能な形であり、しかも「積和交換」操作後は帰納的に消去可能な形にちょうどなってくれるということ。と考えれば何も不思議なところはないよ、と突っ込むべきであったが、議論の始まりをぼ〜っとして聞いてなくて、理解したときにはほぼ収束していた。
- 8.97 式の左辺には ln が必要(あるいは右辺に exp)
- 8.4.6 以降は一応担当外だったけど、"Factor graphs and the sum-product algorithm" という tutorial paper がおもしろかったよ、Turbo code とか LDPC code とか、と付け焼き刃で紹介していたら、以前 LDPC やってたという id:tsubosaka さんにツッコミをもらって、さらに LDPC について簡単ながらちゃんとした紹介をしてもらった。いい感じ。
次回は 3/7(Sun)、９章 EM アルゴリズム＋１０章初め。９章は読んでてワクワクしまくりなので、今から楽しみ。
って、もう定員半分以上埋ってんだけど……
R で K-means クラスタリング
で、最近Ｒ触ってないな〜 K-means 簡単そうだな〜と思ったので、リハビリがてらにさっくり書いてみた。
PRML ９章で使われている Old Faithful のデータセットは PRML サポートサイトにて配布されているのでダウンロードして……なあんてしなくても、R に最初っから入っている。
# Old Faithful dataset を取得して正規化 data("faithful"); xx <- scale(faithful, apply(faithful, 2, mean), apply(faithful, 2, sd));
縦横のスケールが違いすぎて、ユークリッド距離な K-means にこのままでは食わせられないので、正規化してある。
続いて K-means。
> kmeans(xx, 2) K-means clustering with 2 clusters of sizes 98, 174 Cluster means: eruptions waiting 1 -1.2577669 -1.1993566 2 0.7083975 0.6754997 Clustering vector: 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 2 1 2 1 2 1 2 2 1 2 1 2 2 1 2 1 1 2 1 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 2 1 1 2 2 2 2 1 2 2 2 2 2 2 2 2 1 1 2 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 1 2 2 1 2 1 2 2 2 1 2 1 2 2 1 2 1 2 2 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 1 2 2 1 2 1 2 1 2 2 2 1 2 2 1 2 2 1 2 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 1 2 2 2 2 2 2 1 2 2 2 2 1 2 1 2 1 2 1 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 2 2 2 1 2 1 2 1 2 2 1 2 1 2 2 2 1 2 2 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 1 2 1 2 1 2 1 2 2 1 2 2 1 2 1 2 1 2 1 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 2 1 2 1 2 1 2 2 1 2 2 2 1 2 1 2 1 2 2 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 1 2 2 2 2 2 1 2 1 2 1 2 2 2 1 2 1 2 1 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 1 2 2 2 2 2 1 2 2 1 2 2 2 1 2 2 1 2 1 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 2 1 2 2 2 2 2 2 1 2 1 2 2 1 2 1 2 2 1 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 2 2 2 2 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 2 2 2 1 2 1 2 1 1 2 2 1 2 1 2 1 2 2 1 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 2 1 2 1 2 2 2 2 2 2 2 1 2 2 2 1 2 1 1 267 268 269 270 271 272 2 2 1 2 1 2 Within cluster sum of squares by cluster: [1] 24.89206 54.39134 Available components: [1] "cluster" "centers" "withinss" "size"
……すいません、冗談です。
まじめに。
# クラス数 K <- 2; # 中心の初期値(正規乱数) mu <- matrix(rnorm(K*ncol(xx)), nrow=K); # 中心の初期値(本と同じように繰り返し回数が多くなる点を恣意的に指定) mu <- matrix(c(1,-1,-1,1),2)
そして K-means の繰り返し処理。
今回は「式一つだけで K-means の繰り返し処理を書く」というのにこだわってみた。
ノルムを取るところはもうちょっと短く書ける気もしてきたが、もういいかな、と(苦笑
# K-means の１ステップ (mu <- t(sapply(1:K,function(k)colMeans(xx[max.col(-sapply(1:K,function(i)apply(xx,1,function(x)sum((mu[i,]-x)^2))))==k,]))));
実行の様子。
> K <- 2; > mu <- matrix(rnorm(K*ncol(xx)), nrow=K); > (mu <- t(sapply(1:K,function(k)colMeans(xx[max.col(-sapply(1:K,function(i)apply(xx,1,function(x)sum((mu[i,]-x)^2))))==k,])))); eruptions waiting [1,] 0.7400997 1.2204005 [2,] -0.1389707 -0.2291582 > (mu <- t(sapply(1:K,function(k)colMeans(xx[max.col(-sapply(1:K,function(i)apply(xx,1,function(x)sum((mu[i,]-x)^2))))==k,])))); eruptions waiting [1,] 0.7918391 0.7943739 [2,] -0.9040683 -0.9069624 > (mu <- t(sapply(1:K,function(k)colMeans(xx[max.col(-sapply(1:K,function(i)apply(xx,1,function(x)sum((mu[i,]-x)^2))))==k,])))); eruptions waiting [1,] 0.731113 0.7048363 [2,] -1.199593 -1.1564790 > (mu <- t(sapply(1:K,function(k)colMeans(xx[max.col(-sapply(1:K,function(i)apply(xx,1,function(x)sum((mu[i,]-x)^2))))==k,])))); eruptions waiting [1,] 0.7193933 0.6862599 [2,] -1.2373565 -1.1803670 > (mu <- t(sapply(1:K,function(k)colMeans(xx[max.col(-sapply(1:K,function(i)apply(xx,1,function(x)sum((mu[i,]-x)^2))))==k,])))); eruptions waiting [1,] 0.7146233 0.680211 [2,] -1.2487861 -1.188651 > (mu <- t(sapply(1:K,function(k)colMeans(xx[max.col(-sapply(1:K,function(i)apply(xx,1,function(x)sum((mu[i,]-x)^2))))==k,])))); eruptions waiting [1,] 0.7083975 0.6754997 [2,] -1.2577669 -1.1993566
ちゃんと kmeans() と同じ値に収束している。
グラフを描く。きれいにクラスタリングできている。
# 中心 mu を元に分類グラフを描く plot(xx, col=max.col(-sapply(1:K, function(i)apply(xx, 1, function(x)sum((mu[i,]-x)^2)))), xlab=paste(sprintf("%1.3f",t(mu)),collapse=","), ylab=""); points(mu, col = 1:2, pch = 8)
このグラフの点を見る度に間欠泉を思い描き、死ぬまでに一度は Old Faithful 行ってみたいよなあ、と夢見つつも、今のところは Old Faithful のデータを処理するだけで我慢しておくのでした……
*1:「恒等的に１」であるような自明な因子だって許されうる
