ノート/R/確率モデルによるクラスタリング

サイトトップ

ノート
データ処理
├マイニング
├バスケット解析
├バスケット解析をＲで
├図書貸出をＲで
├モデルクラスタリング
└
テキストマイニング
├テキストマイニングTM
├TMとシソーラス
├PubMedをTM
├TMとMeSH
├TMとNLTK
├テキストマイニングとtagger
├医薬品添付文書ＤＢ
├論文の処理1
├剽窃１
├Bing検索
└
分子進化学
├分子進化学
└
Pythonと論文アクセス
├Pythonを使ってみる
├PythonでPubMed
├Pythonで...続き
└
Macメモ
├Script

&ref(): File not found: "favicon_story.txt" at page "ノート/ノート";

自信ないので
├Winのバッチ
│├スクリプト入門
│├バッチファイル
│├SETコマンド
│├自分のメモ
│└
└

メニュー編集

最新の5件

ノート/ノート
ノート/R
訪問者数　5567　　　　　　最終更新　2012-03-04 (日) 13:42:21

Ｒで確率モデルによるクラスタリング（2010-07-21） 　 †

教科書：　新納浩幸：Rで学ぶクラスタ解析（オーム社、ISBN978-4-274-06703-7）の第６章
　　　＞金明哲先生：「Rによるデータサイエンス - データ解析の基礎から最新手法まで」　　＞アマゾンのページ
　　　＞金明哲先生の教材ページ　　＞シリーズのトップページ
　　　＞データマイニング+WEB勉強会のページ

　　　＞Ｒ用パッケージMclustのページ（Washington大 統計学部）　（http://www.stat.washington.edu/mclust/）

　　　＞ＲではないGMM (Gaussian Mixture Model) のパッケージ　〜　mlpack（machine learning package, C++のパッケージ）中のgmmライブラリ
　　　　　Apache Mahoutとそのアルゴリズム
　　　　　Shogun toolbox

理論編（混合分布とそのEMアルゴリズムによる推定） †

• 混合分布
  • 教科書（１）「データマイニングの基礎」（元田浩etal）　オーム社　3.2の6　(p85)
  • 教科書（２）「パターン認識と機械学習」（C.M.ビショップ、元田浩etal訳）　Springer　第９章

• 金田尚久・新居玄武：混合分布問題　その基礎からカーネル降下法まで

• 正田備也：ガウス混合分布のEMアルゴリズムによるパラメータ推定

• 田中洸一：混合分布モデル

• 江口真透：混合分布

• 田中研太郎:混合分布モデルにおける一致推定量の構成

• EMアルゴリズムを用いた混合分布の推定

• 赤穂昭太郎 ECM 法を用いた確率分布の位置, 尺度, 回転パラメータの推定法

• （計算パッケージ）
  • 混合分布に対する EM アルゴリズム (テスト バージョン)
  • OpenCVにあるEMアルゴリズム:EMアルゴリズム（Expectation-Maximization）

• sage/PRML - 混合ガウス分布のEMアルゴリズム

• Tutorial of the EM Algorithm and its Application to Parameter Estimation for Gaussian Mixture and Hidden Markov Models

• Dr. Dowe's page on Mixture Modelling

• EMアルゴリズムの改良、　http://d.hatena.ne.jp/isseing333/20110412/1302619320

• 名古屋大能登原

• 北大井上

• 高知大本田

• 茨城大新納

• EMアルゴリズムの教科書

• 機械学習

• 多変量混合正規分布のEMアルゴリズム・サンプルプログラム　その１

Ｒによる解析５ （山の具合を描いてみた）　（2010/07/21　17:30） †

山の形に見えるのが結構難しい。
３次元プロットすると、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）しか自動クラスタ化をしていないので、それで比較する。

左側の図から山を読み取ると、

• X=36付近に高さ80程度のピーク、Xを1-1024のスケールに直すと16倍してX=576。これは右図の緑に対応するだろう。
• X=38付近に高さ80程度のピーク、Xを1-1024のスケールに直すと16倍してX=608。これは右図の青に対応するだろう。
• X=42〜49付近に高さ30程度のピーク、Xを1-1024のスケールに直すと16倍してX=672〜784。これは右図の赤に対応するだろう。

Ｒによる解析４ （FL1もクラスタ化してみた）　（2010/07/16　18:30） †

FSC-FL1は、自動クラスタ化で、クラスタ数は６になった。

　　　

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

Ｒによる解析３（点数を増やしてみると）　(2010/07/16） †

[A] クラスタ数自動決定(=8)の場合 †

計算時間を余分にかけて、測定点数を増やしてみる。同一実験データ（100113 24h 40C +20-1_nonlimit.csv）で、10,000点から20,000点へ。クラスタ数自動決定をすると８になった。

　　　

[B] クラスタ数固定（=6）の場合 †

一方で、クラスタ数を10,000点時と同じ６つに固定すると、次のようなクラスタを作る。

　　　

[A] クラスタ数自動決定の場合の出力データ †

20,000点時、クラスタ数＝８（自動決定）の出力データ

$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

[B] クラスタ数を6に固定した場合の出力データ †

クラスタ数を10,000点時と同じ６つに固定したときのパラメタ出力。

$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

Ｒによる解析２（雑音除去の効果と、Mclustでの表示法の改善）　（2010/07/14） †

雑音除去がかなり効果を挙げることがわかった。具体的には、(1)プロットチャートの左端の線（Xの値が0）と、(2)全体に散らばっている≪おそらくランダムノイズと思われる点≫、の２種類を除去すると、余分なヤマが出なくなる。

例：100113 24h 40C +20-1_nonlimit.csvを入力とする場合

• 全体で20,000個のデータがあるが、時間短縮のため最初の10,000個のデータのみ利用した
• それに対して、上記の(1),(2)での除去を行うと結果は2,592個のデータが残った。かなり少なくなったと言えるだろう。
• その2,592個のデータに対して、Mclustを用いて（正規分布として）ヤマを推定した。
• 結果は

  最適のモデルはVVV、ヤマの数は６つ。
  
  それぞれの分布（ヤマ）の割合
  $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)

のような関数を使うことが出来る。結果は、

色の「凡例」が表示されていないのでわかりづらいが、
　　グループ１＝青、２＝赤、３＝黄緑、４＝空色、、５＝ピンク、６＝濃い緑
に対応する。

真ん中の赤（グループ２、データ点数538）と青（グループ１、データ点数942）がメイン、 黄緑（グループ３、データ点数622）が右上に新しく伸びた別の山、 右上の濃い緑（グループ６、データ点数322）が定点用ビーズ、 更に右上の空色（グループ４、データ点数92）と左下のピンク（グループ５、データ点数76） はたぶんノイズと思われる。

このように、かなりきれいな形で分布と測定点の帰属が推定できると思われる。他のデータについてもテストしてみたい。

ゼミで使った資料　　（2010/07/12追記） †

Ｒによる解析（新納の6.4節）（2010/06/29） †

手順

• パッケージmclustをインストール
  　　パッケージ　→　パッケージのインストール　→　mclustを選択
• プログラム例は

  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データ＝からはじめの４値を抜き出したもの。

  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
  ...

フローサイトメトリー（大腸菌）データでのモデルベースクラスタリング †

実験１　〜　測定データ　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

となっている。大きいクラスタとしてクラスタ２（551,578）、クラスタ６（609,581）、クラスタ５の(722,701)を見ると、分布図（x-y座標とも、数値を1/4にしている＝1024値を256等分、また前回の処理の都合上FSCの100以下、SSCの50以下を切り取ってしまっていることに注意、クラスタ推定は切り取り前のデータを利用）上の点(138,145)、 （152,145）、 （181,175）のピーク（前者は赤い点２つ、最後の点は右上に伸びた塊の頂上）と対応付けられる。

この解析では、それぞれの分布に属する点の数を、モデルの形から推定できていることであり、具体的には上記のclassification tableに相当する。6000点の中からクラスタ５（右上に伸びているクラスタ）に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