関数クラスタリングのMFclust関数の内部構造

関数クラスタリングを行うMFclust関数の内部構造を調査しました（MDFAパッケージ）。

この手法は遺伝子系の論文Nucleic Acids Researchに掲載された手法です。
http://genemerge.cbcb.umd.edu/online/SSC.pdf

ある遺伝子の経時変化を関数クラスタリングで分類して、クラスタ毎に曲線を当てはめるというもの。

MFclustの内容は次のようになっています。

• 全個体、全時系列データに対してk-means法を当てはめる
  • 時系列のパターンによってクラスタ分けされる

• それぞれのクラスタに対して時間軸をxとしたSmoothing Spline ANOVA Models（スプラインANOVA）を当てはめる
  • 関数ANOVA（functional ANOVA）の一種
  • このとき個体（個人）を変量効果に入れる
  • 次の図のようなイメージ

• パラメータ推定をEMアルゴリズムで行う

• 最初のk-means法のクラスタ数をいくつか試し、それぞれのクラスタ数でのBICを計算する
  • BICが最小のクラスタ数が最適なクラスタになる

• 最適なクラスタ数での、それぞれのクラスタの時系列関数が求まる

このように、今までブログで紹介してきた変量効果モデルやEMアルゴリスムなどがミックスされています。

まるでこの論文のために今まで紹介してきたような感覚にさえ陥りますorz

それはさておき、、MFclust関数は内部で他の関数を呼び出しています。

1つの関数で全てを描くとメンテが大変なので分けているのですが、それらの関数の関係を次に紹介します。

MFclust
┣MFclust.compute （minG〜maxGの間で繰り返す）
┃┃#初期値の計算
┃┣ kmeans
┃┣compute.M.all （kmeansのクラスの数だけ繰り返す）
┃┃┗compute.center
┃┃ ┗ssanova （Smoothing Spline ANOVA Models、変量効果も指定）
┃┃#初期値を使ってEMアルゴリズム
┃┗em.clust （iter.maxの数だけ繰り返す）
┃ ┃#Eステップ
┃ ┣Estep.tk
┃ ┃#Mステップ
┃ ┣compute.weight
┃ ┣compute.M.pk
┃ ┣compute.M.all （weightを指定する）
┃ ┃#BICの計算
┃ ┗em.bic
┃#BICが最小のクラスター数を選択する
┗min

このように、多くの関数で成り立っています。

大まかに推定の手順を書き下します。

1. k-meansでクラスタ分類
2. Mステップを行って初期値を求める
   • ssanova関数に変量効果を入れてパラメータを推定
3. Eステップを行って期待値を計算する
   • Mステップでは「重み」として利用する
4. Mステップで目的関数を最大化する
   • 重み付きのssanovaを実行してパラメータを推定
5. EステップとMステップを繰り返す
6. 尤度が最大のときの値を最適解として、そのときのBICを求める
7. kを動かしていくつかBICを求め、最小のときのクラスタ数を選択する

EMアルゴリズムと大層な事を言っても、やってるのはssanova関数の当てはめ（Mステップ）と重みの計算（Eステップ）。

こうやって関数の中身を見ると実装のイメージが付きやすいですね。

論文中ではこのように説明されてます。

しかし自分でこれを書くのはやっぱり大変です。

これくらいの方法を考え、実装できると文句なしに博士号だと思います！