![]() |
ノート/R/確率モデルによるクラスタリングhttps://pepper.is.sci.toho-u.ac.jp:443/pepper/index.php?%A5%CE%A1%BC%A5%C8%2FR%2F%B3%CE%CE%A8%A5%E2%A5%C7%A5%EB%A4%CB%A4%E8%A4%EB%A5%AF%A5%E9%A5%B9%A5%BF%A5%EA%A5%F3%A5%B0 |
![]() |
ノート/ノート
ノート/R
訪問者数 5755 最終更新 2012-03-04 (日) 13:42:21
教科書: 新納浩幸:Rで学ぶクラスタ解析(オーム社、ISBN978-4-274-06703-7)の第6章
>金明哲先生:「Rによるデータサイエンス - データ解析の基礎から最新手法まで」 >アマゾンのページ
>金明哲先生の教材ページ >シリーズのトップページ
>データマイニング+WEB勉強会のページ
>R用パッケージMclustのページ(Washington大 統計学部) (http://www.stat.washington.edu/mclust/)
>RではないGMM (Gaussian Mixture Model) のパッケージ 〜 mlpack(machine learning package, C++のパッケージ)中のgmmライブラリ
Apache Mahoutとそのアルゴリズム
Shogun toolbox
山の形に見えるのが結構難しい。
3次元プロットすると、X-Yの解像度が高いと、細い縦の線になってしまう。つまり表面が平滑化されていない。そこで、X-Yの解像度を256, 128, 64で比べてみる。
[A] 少ないデータの方(FSC-SSC_100113_24h_40C_+20-1)
[B] 多いデータの方(FSC-SSC_100113_24h_40C_+9-1)
[C]自動クラスタ化の結果との比較
上記を踏まえて、X-Y解像度を64に設定した図で考えることにする。
少ないデータの方(FSC-SSC_100113_24h_40C_+20-1)しか自動クラスタ化をしていないので、それで比較する。
左側の図から山を読み取ると、
FSC-FL1は、自動クラスタ化で、クラスタ数は6になった。
総データ数 | FSC-FL1 | FSC-SSC | |||||||||||
クラスタ番号 | 1 | 2 | 3 | 4 | 5 | 6 | 1 | 2 | 3 | 4 | 5 | 6 | |
(色別) | 青 | 赤 | 黄緑 | 空色 | ピンク | 濃緑 | 青 | 赤 | 黄緑 | 空色 | ピンク | 濃緑 | |
両図のクラスタ間の対応は未決定なので、要注意 | |||||||||||||
クラスタの大きさ(割合) | 0.23 | 0.27 | 0.30 | 0.42 | 0.03 | 0.13 | 0.37 | 0.29 | 0.21 | 0.02 | 0.02 | 0.08 | |
同上 (点数) | 752 | 876 | 973 | 135 | 92 | 417 | 1694 | 1308 | 951 | 102 | 110 | 398 | |
クラスタの平均 x | 712 | 617 | 567 | 998 | 408 | 904 | 602 | 690 | 545 | 999 | 376 | 904 | |
同上 y | 631 | 531 | 514 | 978 | 289 | 981 | 586 | 666 | 560 | 1012 | 253 | 941 | |
クラスタの分散 x-x | 2129 | 951 | 2375 | 164 | 9637 | 18 | 947 | 3150 | 2791 | 190 | 6395 | 14 | |
同上 y-x | 1654 | 983 | 1046 | -1 | -274 | -1 | 321 | 2193 | 937 | 106 | -1959 | 4 | |
同上 y-y | 1841 | 2411 | 1598 | 1 | 15453 | 2 | 385 | 1983 | 990 | 63 | 5303 | 7 |
FSC-FL1の出力データの詳細
$pro [1] 0.23167330 0.26984454 0.29997641 0.04160247 0.02839789 0.12850539 $mean [,1] [,2] [,3] [,4] [,5] [,6] [1,] 712.3932 616.5135 567.3878 998.4148 408.2780 903.6547 [2,] 631.3966 531.1113 514.4008 978.0889 288.7181 981.1727 $variance $variance$modelName [1] "VVV" $variance$d [1] 2 $variance$G [1] 6 $variance$sigma 1 [,1] [,2] [1,] 2129.354 1654.443 [2,] 1654.443 1840.998 2 [,1] [,2] [1,] 950.7647 982.5923 [2,] 982.5923 2411.1313 3 [,1] [,2] [1,] 2374.617 1046.110 [2,] 1046.110 1597.614 4 [,1] [,2] [1,] 163.813114 -1.2368724 [2,] -1.236872 0.4661728 5 [,1] [,2] [1,] 9637.2769 -273.5086 [2,] -273.5086 15453.0124 6 [,1] [,2] [1,] 17.5594086 -0.6310232 [2,] -0.6310232 2.3634733 [1] "Total data counts: 6490" [1] "Cluster sizes:" [1] 1504 1751 1947 270 184 834 [,1] [,2] [,3] [,4] [,5] [,6] [1,] 712.3932 616.5135 567.3878 998.4148 408.2780 903.6547 [2,] 631.3966 531.1113 514.4008 978.0889 288.7181 981.1727
計算時間を余分にかけて、測定点数を増やしてみる。同一実験データ(100113 24h 40C +20-1_nonlimit.csv)で、10,000点から20,000点へ。クラスタ数自動決定をすると8になった。
総データ数 | 20,000 | 10,000 | |||||||||||||
クラスタ番号 | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 1 | 2 | 3 | 4 | 5 | 6 | |
(色別) | 青 | 赤 | 黄緑 | 空色 | ピンク | 濃緑 | 紫 | 黄土 | 青 | 赤 | 黄緑 | 空色 | ピンク | 濃緑 | |
クラスタの大きさ(割合) | 0.11 | 0.31 | 0.27 | 0.11 | 0.02 | 0.02 | 0.09 | 0.07 | 0.36 | 0.21 | 0.24 | 0.04 | 0.03 | 0.12 | |
同上 (点数) | 1028 | 2823 | 2421 | 1014 | 204 | 217 | 796 | 621 | 942 | 538 | 622 | 92 | 76 | 322 | |
クラスタの平均 x | 619 | 586 | 659 | 524 | 999 | 377 | 903 | 748 | 608 | 547 | 699 | 1000 | 357 | 904 | |
同上 y | 580 | 585 | 643 | 546 | 1012 | 249 | 941 | 700 | 589 | 568 | 674 | 1013 | 255 | 941 | |
クラスタの分散 x-x | 526 | 1234 | 2978 | 3026 | 190 | 6398 | 14 | 966 | 582 | 3426 | 2525 | 189 | 7504 | 13 | |
同上 y-x | 468 | 361 | 2377 | 775 | 106 | -1724 | 4 | 899 | 250 | 1041 | 1743 | 106 | -1003 | 4 | |
同上 y-y | 537 | 345 | 2438 | 875 | 63 | 4444 | 7 | 989 | 353 | 931 | 1677 | 64 | 2800 | 6 |
一方で、クラスタ数を10,000点時と同じ6つに固定すると、次のようなクラスタを作る。
総データ数 | 20,000 | 10,000 | |||||||||||
クラスタ番号 | 1 | 2 | 3 | 4 | 5 | 6 | 1 | 2 | 3 | 4 | 5 | 6 | |
(色別) | 青 | 赤 | 黄緑 | 空色 | ピンク | 濃緑 | 青 | 赤 | 黄緑 | 空色 | ピンク | 濃緑 | |
クラスタ2(赤)と3(黄緑)が逆転しているので注意 | |||||||||||||
クラスタの大きさ(割合) | 0.37 | 0.29 | 0.21 | 0.02 | 0.02 | 0.08 | 0.36 | 0.21 | 0.24 | 0.04 | 0.03 | 0.12 | |
同上 (点数) | 1694 | 1308 | 951 | 102 | 110 | 398 | 942 | 538 | 622 | 92 | 76 | 322 | |
クラスタの平均 x | 602 | 690 | 545 | 999 | 376 | 904 | 608 | 547 | 699 | 1000 | 357 | 904 | |
同上 y | 586 | 666 | 560 | 1012 | 253 | 941 | 589 | 568 | 674 | 1013 | 255 | 941 | |
クラスタの分散 x-x | 947 | 3150 | 2791 | 190 | 6395 | 14 | 582 | 3426 | 2525 | 189 | 7504 | 13 | |
同上 y-x | 321 | 2193 | 937 | 106 | -1959 | 4 | 250 | 1041 | 1743 | 106 | -1003 | 4 | |
同上 y-y | 385 | 1983 | 990 | 63 | 5303 | 7 | 353 | 931 | 1677 | 64 | 2800 | 6 |
20,000点時、クラスタ数=8(自動決定)の出力データ
$pro [1] 0.11269930 0.30937046 0.26537431 0.11109090 0.02235861 0.02375630 [7] 0.08724244 0.06810768 $mean [,1] [,2] [,3] [,4] [,5] [,6] [,7] [1,] 619.0596 585.5307 658.6493 523.5670 998.6863 376.6446 903.6633 [2,] 579.8497 584.8812 643.2790 546.4888 1011.5196 248.8440 940.8317 [,8] [1,] 748.1149 [2,] 700.2097 $variance$modelName [1] "VVV" $variance$d [1] 2 $variance$G [1] 8 $variance$sigma 1 [,1] [,2] [1,] 526.1028 467.6122 [2,] 467.6122 537.2284 2 [,1] [,2] [1,] 1234.2949 361.0131 [2,] 361.0131 345.3007 3 [,1] [,2] [1,] 2977.924 2376.559 [2,] 2376.559 2437.739 4 [,1] [,2] [1,] 3026.429 774.5300 [2,] 774.530 874.8733 5 [,1] [,2] [1,] 190.1761 105.9473 [2,] 105.9473 63.2104 6 [,1] [,2] [1,] 6397.935 -1723.837 [2,] -1723.837 4444.096 7 [,1] [,2] [1,] 14.062524 3.935784 [2,] 3.935784 6.748043 8 [,1] [,2] [1,] 965.6994 899.0406 [2,] 899.0406 988.8561 [1] "Total data counts: 9124" [1] "Cluster sizes:" [1] 1028 2823 2421 1014 204 217 796 621 [,1] [,2] [,3] [,4] [,5] [,6] [,7] [1,] 619.0596 585.5307 658.6493 523.5670 998.6863 376.6446 903.6633 [2,] 579.8497 584.8812 643.2790 546.4888 1011.5196 248.8440 940.8317 [,8] [1,] 748.1149 [2,] 700.2097
クラスタ数を10,000点時と同じ6つに固定したときのパラメタ出力。
$pro [1] 0.37127726 0.28670306 0.20830639 0.02235861 0.02411224 0.08724244 $mean [,1] [,2] [,3] [,4] [,5] [,6] [1,] 601.7230 689.6276 545.2734 998.6863 375.6309 903.6633 [2,] 585.6392 666.4851 560.2669 1011.5196 252.5498 940.8317 $variance $variance$modelName [1] "VVV" $variance$d [1] 2 $variance$G [1] 6 $variance$sigma 1 [,1] [,2] [1,] 947.0633 321.2171 [2,] 321.2171 385.4888 2 [,1] [,2] [1,] 3149.918 2193.402 [2,] 2193.402 1983.022 3 [,1] [,2] [1,] 2790.827 936.6290 [2,] 936.629 990.1658 4 [,1] [,2] [1,] 190.1761 105.9473 [2,] 105.9473 63.2104 5 [,1] [,2] [1,] 6394.872 -1959.228 [2,] -1959.228 5302.817 6 [,1] [,2] [1,] 14.062524 3.935784 [2,] 3.935784 6.748043 [1] "Total data counts: 9124" [1] "Cluster sizes:"
[1] 1694 1308 951 102 110 398 [,1] [,2] [,3] [,4] [,5] [,6] [1,] 601.7230 689.6276 545.2734 998.6863 375.6309 903.6633 [2,] 585.6392 666.4851 560.2669 1011.5196 252.5498 940.8317
雑音除去がかなり効果を挙げることがわかった。具体的には、(1)プロットチャートの左端の線(Xの値が0)と、(2)全体に散らばっている≪おそらくランダムノイズと思われる点≫、の2種類を除去すると、余分なヤマが出なくなる。
例:100113 24h 40C +20-1_nonlimit.csvを入力とする場合
最適のモデルはVVV、ヤマの数は6つ。 それぞれの分布(ヤマ)の割合 $pro [1] 0.36346510 0.20738490 0.24009412 0.03549383 0.02933365 0.12422840 それぞれの分布(ヤマ)のデータ個数 (2,592×$pro) [1] 942 538 622 92 76 322 それぞれの分布(ヤマ)の中心(平均値) [,1] [,2] [,3] [,4] [,5] [,6] [1,] 608.3974 547.3208 699.312 999.8478 356.7617 903.913 [2,] 589.3498 567.8984 673.731 1012.5870 254.5305 941.149 それぞれの分布(ヤマ)の分散(sigma) 1 2 3 [,1] [,2] [,1] [,2] [,1] [,2] [1,] 581.7725 250.0036 [1,] 3425.921 1040.665 [1,] 2524.561 1742.984 [2,] 250.0036 353.2514 [2,] 1040.665 931.292 [2,] 1742.984 1677.382 4 5 6 [,1] [,2] [,1] [,2] [,1] [,2] [1,] 188.6073 106.19802 [1,] 7504.097 -1002.931 [1,] 12.849581 4.329733 [2,] 106.1980 63.98157 [2,] -1002.931 2800.267 [2,] 4.329733 6.114425
また、表示法として、
mclust2Dplot(data=z1, what="classification", identify=TRUE, parameters=x$paramters, z=x$z, scale=TRUE)
のような関数を使うことが出来る。結果は、
色の「凡例」が表示されていないのでわかりづらいが、
グループ1=青、2=赤、3=黄緑、4=空色、、5=ピンク、6=濃い緑
に対応する。
真ん中の赤(グループ2、データ点数538)と青(グループ1、データ点数942)がメイン、 黄緑(グループ3、データ点数622)が右上に新しく伸びた別の山、 右上の濃い緑(グループ6、データ点数322)が定点用ビーズ、 更に右上の空色(グループ4、データ点数92)と左下のピンク(グループ5、データ点数76) はたぶんノイズと思われる。
このように、かなりきれいな形で分布と測定点の帰属が推定できると思われる。他のデータについてもテストしてみたい。
手順
library(mclust) x <- matrix(scan("myiris.dat"), ncol=4, byrow=TRUE) mc <- Mclust(x, G=3,modelNames="EII") ans <- mc$classification pa <- mc$parameters pa$pro pa$mean plot(Mclust(x, G=1:20,modelNames="VEI")) # modelNames="EEI","VII", etcここで、myiris.datはR中のサンプルデータiris= (花びらの長さ・幅・がく片の長さ・幅・irisの種類)が150データ=からはじめの4値を抜き出したもの。
5.1 3.5 1.4 0.2 4.9 3.0 1.4 0.2 4.7 3.2 1.3 0.2 4.6 3.1 1.5 0.2 5.0 3.6 1.4 0.2 ...
実験1 〜 測定データ FSC-SSCをそのまま入力してモデルベースのクラスタリングをしてみる。
データはほとんどそのまま使える。
install.packages("mclust") library(mclust) # Read the data ##u <- matrix(scan("100113 24h 40C +9-1_nonlimit.csv", sep=","), ncol=8, byrow=TRUE) u <- matrix(scan("100113 24h 40C +20-1_nonlimit.csv", sep=","), ncol=8, byrow=TRUE) # ncol should be 8 because lines end with comma (with one more hidden column.) u <- u[,1:7] # erase the hidden=empty column
mclustの使い方については、いろいろ書いてあるところもあるが、全体を見るにはオリジナルのドキュメントも良い。--> tr504.pdf
# The document starts with BIC <- mclustBIC(u[1:6000,1:2])
uのうちの1:2 (FSCとSSC)を使ってみる。また、全データ(10000点)を使うとメモリ不足になるので、最初の4000点を使ってみる。Windowsで6000点ぐらいまでは動いた。なお、メモリスペースについては
help(memory.size)
で説明が出てくるが、memory.size(NA) と打つと現在のメモリサイズ(Windows XPの場合デフォルトで2GB、最大で3GB、Win7では4GBらしい)が表示され、
memory.limit(3071)
とすると上限を3GBに設定することが出来る。いずれにせよ、この程度では、全ポイントを含めると、メモリが足りなくなる。
結果で得られたBICを印刷すると以下の通りで、プロットすると図の通り。
BIC: EII VII EEI VEI EVI VVI EEE 1 -162779.0 -162779.0 -162262.9 -162262.9 -162262.9 -162262.9 -155018.1 2 -156268.7 -155705.2 -155662.7 -155286.2 -155283.2 -155294.3 -153774.0 3 -154328.1 -153104.1 -154337.0 -153033.5 -153878.3 NA -149468.2 4 -151256.3 -149930.7 -151269.8 -149849.3 -150173.3 NA -148455.5 5 -149347.0 -147278.0 -149355.8 -147088.4 -147818.7 NA -148300.4 6 -148476.1 -145945.0 -148481.6 -145875.4 -146298.5 NA -147772.1 7 -148497.6 -144866.6 -148494.2 -144805.9 -145657.5 NA -147800.0 8 -148454.1 -144143.6 -148467.9 -144108.7 -145461.2 NA -147827.8 9 -146690.5 -143815.6 -146592.5 -143822.5 -144603.3 NA -147842.7 EEV VEV VVV 1 -155018.1 -155018.1 -155018.1 2 -150277.4 -148559.5 -148223.5 3 -147289.1 -145205.7 NA 4 -145856.6 -144873.9 NA 5 -144800.7 -143775.9 NA 6 -144804.5 -142855.7 NA 7 -144002.2 -142175.9 NA 8 -143779.0 -142215.5 NA 9 -143724.4 -141796.0 NA
要するに、方法はVEVが最もよくて、クラスター数は9ということである。
出力した結果のオリジナル表記は下記の通りになっていた。
> summary(BIC, data=u[1:6000,1:2],modelNames="VEV") classification table: 1 2 3 4 5 6 7 8 9 117 1231 639 431 1013 1400 443 591 135 best BIC values: VEV,9 VEV,7 VEV,8 -141796.0 -142175.9 -142215.5
となっていた。再度
x <- Mclust(u[1:6000,1:2], G=9, modelNames="VEV")
を処理して、推定したパラメタを見てみると、まずそれぞれのクラスタの大きさαiは
$pro [1] 0.06943546 0.20443535 0.09925967 0.09049414 0.13425212 0.17498290 [7] 0.10774480 0.09694017 0.02245540
で、i=2, 3, 5, 6, 8などが大きい。次にそれぞれのクラスタの平均値と分散は、
$mean [,1] [,2] [,3] [,4] [,5] [,6] [,7] [1,] 672.5915 550.8010 350.8723 348.1141 721.9350 609.336 477.4751 [2,] 652.4011 578.4403 238.2852 560.2122 701.1421 581.061 466.0265 [,8] [,9] [1,] 0.3521280 903.9811 [2,] 252.0523710 941.5063 $variance$sigma , , 1 [,1] [,2] [1,] 4510.760 2559.932 [2,] 2559.932 1932.424 , , 2 [,1] [,2] [1,] 5840.583 2278.554 [2,] 2278.554 1350.294 , , 3 [,1] [,2] [1,] 18864.03112 -47.12685 [2,] -47.12685 1101.31457 , , 4 [,1] [,2] [1,] 120829.0 57238.70 [2,] 57238.7 37994.28 , , 5 [,1] [,2] [1,] 3609.556 2464.675 [2,] 2464.675 2164.224 , , 6 [,1] [,2] [1,] 2413.461 1707.167 [2,] 1707.167 1546.101 , , 7 [,1] [,2] [1,] 15796.04 14299.09 [2,] 14299.09 16354.06 , , 8 [,1] [,2] [1,] 84.9818781 -0.9859035 [2,] -0.9859035 1455.7868607 , , 9 [,1] [,2] [1,] 47.08964 19.93518 [2,] 19.93518 12.34935
となっている。大きいクラスタとしてクラスタ2(551,578)、クラスタ6(609,581)、クラスタ5の(722,701)を見ると、分布図(x-y座標とも、数値を1/4にしている=1024値を256等分、また前回の処理の都合上FSCの100以下、SSCの50以下を切り取ってしまっていることに注意、クラスタ推定は切り取り前のデータを利用)上の点(138,145)、 (152,145)、 (181,175)のピーク(前者は赤い点2つ、最後の点は右上に伸びた塊の頂上)と対応付けられる。
この解析では、それぞれの分布に属する点の数を、モデルの形から推定できていることであり、具体的には上記のclassification tableに相当する。6000点の中からクラスタ5(右上に伸びているクラスタ)に1013点が属していることが分かる。
その他の出力は
$variance$scale [1] 1470.86706 1641.55235 4557.74192 36256.59046 1318.05373 903.89915 [7] 7339.30971 351.73076 13.56889 $variance$shape [1] 4.1389259 0.2416086 $variance$orientation , , 1 [,1] [,2] [1,] 0.8514047 -0.5245093 [2,] 0.5245093 0.8514047 , , 2 [,1] [,2] [1,] 0.9224601 0.3860924 [2,] 0.3860924 -0.9224601 , , 3 [,1] [,2] [1,] -0.999996481 -0.002653105 [2,] 0.002653105 -0.999996481 , , 4 [,1] [,2] [1,] -0.8905668 -0.4548525 [2,] -0.4548525 0.8905668 , , 5 [,1] [,2] [1,] 0.8004262 0.5994314 [2,] 0.5994314 -0.8004262 , , 6 [,1] [,2] [1,] 0.7893714 0.6139159 [2,] 0.6139159 -0.7893714 , , 7 [,1] [,2] [1,] 0.7001755 0.7139708 [2,] 0.7139708 -0.7001755 , , 8 [,1] [,2] [1,] -0.0007192144 -0.9999997414 [2,] 0.9999997414 -0.0007192144 , , 9 [,1] [,2] [1,] -0.9102023 -0.4141639 [2,] -0.4141639 0.9102023