- 追加された行はこの色です。
- 削除された行はこの色です。
#freeze
*ベータ混合分布モデル [#e48dc134]
#contents
**ベータ分布とは／単独のベータ分布推定 [#dc72364f]
-[[MATLAB Note/統計/確率分布/ベータ分布（Beta distribution）:http://shower.human.waseda.ac.jp/~m-kouki/pukiwiki_public/115.html#ff6a4cd6]] を参照。
-[[MATLAB Note/統計/確率分布/ベータ分布（Beta distribution）:http://speechresearch.fiw-web.net/115.html#ff6a4cd6]] を参照。
**複数のベータ分布（ベータ混合分布）の推定 [#pcdc9025]
***目的 [#w83228f3]
-正規分布ではない混合分布のクラスタ数推定を行いたい場合。
--クラスタ数を指定するならGMMでもある程度ロバストだが、BICの結果はデータの正規性の影響を強く受ける。
***どのようなとき使えそうか（MATLABによる予備的検討） [#m910c4b8]
-単独のベータ分布ではフィッティングできないとき
--単独ベータ分布でも両端に山が来るような形状をとれるが、明確な二峰性データのフィッティングは無理。
--例えば &ref(beta_bi.txt); ((これは [[Houseman, Christensen, Yeh, Marsit, Karagas, Wrensch, Nelson, Wiemels, Zheng, Wiencke, & Kelsey, BMC Bioinformatics, 9(365), 2008.:http://www.biomedcentral.com/1471-2105/9/365/]] の caseII で作ったテストデータ（[[Additional files:http://www.biomedcentral.com/1471-2105/9/365/additional]]で公開されているもの）の1列目。))をMATLABでフィットしてみると...
#geshi(matlab){{
data = csvread('beta_bi.txt');
linescale = [0:0.01:1.2];
[phat, pci] = betafit(data) % パラメータ推定
pdfdata = betapdf(linescale, phat(1), phat(2)); % 確率密度関数生成
[histdata, position] = hist(data);
histdata = histdata / (max(histdata) / max(pdfdata));
figure, bar(position, histdata); xlim([0,1.2]); % プロット
hold on
plot(linescale, pdfdata, 'r','LineWidth',2); xlim([0,1.2]); title('betafit');
hold off
}}
#ref(beta_bi.png,,40%);
--このデータは２つのベータ分布の混合（混合ベータ分布）であるので、単独のベータ分布でうまくあてはめることはできなかった。
-ガウス混合分布モデルの当てはまりが悪い（もしくはクラスタ数をうまく推定できない）とき
--上のデータをMATLABのGMMにかけて、AICを求めて最適な混合数も求めてみます。MATLABのstatistics toolboxが必要です。
#geshi(matlab){{
data = csvread('beta_bi.txt');
options = statset('Display','final'); % STATS構造体の表示設定
linescale = [0:0.01:1.2];
for k = 1:4
obj = gmdistribution.fit(data,k,'Options',options); % カテゴリ数k
MU = obj.mu % 平均
SIGMA = obj.Sigma % 共分散
P = obj.PComponents % 混合比
AIC= obj.AIC; % AIC
% 推定したGMMパラメータから混合正規分布を作る
obj_plot = gmdistribution(MU, SIGMA, P);
pdfdata = pdf(obj_plot, linescale'); % 確率密度関数生成
% プロット
subplot(2,2,k);
[histdata, position] = hist(data);
histdata = histdata / (max(histdata) / max(pdfdata)); % 縦軸を揃える
bar(position, histdata); xlim([0,1.2]); % データのヒストグラム
hold on
plot(linescale, pdfdata, 'r','LineWidth',2);
xlim([0,1.2]);
title(strcat('cluster:', int2str(k), ' AIC:', int2str(AIC)));
hold off
end
}}
#ref(beta_bi_GMM.png,,60%);
--AICが最小になるのはクラスタ数が3のときで、正しくクラスタ数を推定できませんでした。また混合正規分布でフィッティングするので正確なデータの分布の形状を再現できていません。
***ベータ混合分布モデルのRスクリプト((CRANに登録されている [[RPMM: Recursively Partitioned Mixture Model:http://cran.r-project.org/web/packages/RPMM/index.html]] も検討したのですが（[[Houseman, et al., 2008.:http://www.biomedcentral.com/1471-2105/9/365/]] の [[Additional file 2:http://www.biomedcentral.com/1471-2105/9/365/additional]] にスクリプトの実行例）、RPMMは超多混合クラスタを分類するために階層型クラスタリングを導入する、というもので、DNAデータの解析に特化しているようでうまくいきませんでした（'''[[検討の跡:http://shower.human.waseda.ac.jp/~m-kouki/pass/E_Andres_Houseman_2008.zip]]'''）)) ((CRANに登録されている [[mixtools: Tools for analyzing finite mixture models:http://cran.r-project.org/web/packages/mixtools/]] も検討したのですがベータ混合分布の推定ができるのかどうかはよくわかっていません（[[検討の過程:http://shower.human.waseda.ac.jp/~m-kouki/mixtools(unsolved).txt]]）)) [#w45470b7]
***ベータ混合分布モデルのRスクリプト((CRANに登録されている [[RPMM: Recursively Partitioned Mixture Model:http://cran.r-project.org/web/packages/RPMM/index.html]] も検討したのですが（[[Houseman, et al., 2008.:http://www.biomedcentral.com/1471-2105/9/365/]] の [[Additional file 2:http://www.biomedcentral.com/1471-2105/9/365/additional]] にスクリプトの実行例）、RPMMは超多混合クラスタを分類するために階層型クラスタリングを導入する、というもので、DNAデータの解析に特化しているようでうまくいきませんでした（'''[[検討の跡:http://shower.human.waseda.ac.jp/~m-kouki/pass/E_Andres_Houseman_2008.zip]]'''（リンク切れ））)) ((CRANに登録されている [[mixtools: Tools for analyzing finite mixture models:http://cran.r-project.org/web/packages/mixtools/]] も検討したのですがベータ混合分布の推定ができるのかどうかはよくわかっていません（[[検討の過程:http://shower.human.waseda.ac.jp/~m-kouki/mixtools(unsolved).txt]]（リンク切れ））)) [#w45470b7]
-[[Yuan Ji 先生:http://odin.mdacc.tmc.edu/~ylji/]] が ''[[betamix_2.R:http://odin.mdacc.tmc.edu/~ylji/betamix_2.R]]'' を公開なさっています。(([[RPMM: Recursively Partitioned Mixture Model:http://cran.r-project.org/web/packages/RPMM/index.html]] の論文（[[Houseman, et al., 2008.:http://www.biomedcentral.com/1471-2105/9/365/]]）で先行研究として引用されています。この論文によれば、どちらも結果の信頼度は同じで、RPMMの方が計算の効率(?)が良いとのこと。))
--そのままでは実行できないのでファイルパスやテストデータなどを一部書き直す必要があります。
---'''[[動くように変更＆日本語注釈:http://shower.human.waseda.ac.jp/~m-kouki/pass/BetaMixtureModel_R.zip]]'''（限定公開です） → '''[[環境研Sさんによる改良版:http://shower.human.waseda.ac.jp/~m-kouki/pass/BetaMixtureModel_R_s.zip]]'''、ありがとうございました！
---'''[[動くように変更＆日本語注釈:http://shower.human.waseda.ac.jp/~m-kouki/pass/BetaMixtureModel_R.zip]]'''（リンク切れ） → '''[[環境研Sさんによる改良版:http://shower.human.waseda.ac.jp/~m-kouki/pass/BetaMixtureModel_R_s.zip]]'''、ありがとうございました！（申し訳ございません、現在リンク切れです）
--[[Ji Y, Wu C, Liu P, Wang J, Coombes KR. Applications of beta-mixture models in bioinformatics. Bioinformatics 21(9):2118-22, 5/2005.:http://bioinformatics.oxfordjournals.org/content/21/9/2118.full.pdf]]
---評価はICL-BIC (Biernacki et al., 1998)
AICやBICはベータ混合分布モデルのクラスタ数の推定には適切でなく、混合数を過大評価しがちである、というのが本論文の新規な指摘点。以下は文献より引用。
One key issue in the finite mixture model approach is to decide the number of components in the model. This has been extensively studied for the normal-mixture model. The commonly used Akaike information criterion (AIC) Akaike (1973) and Bayesian information criterion (BIC) (Schwarz, 1978) have been shown to be adequate for deciding the number of components in the normal-mixture model (Leroux, 1992; Roeder and Wasserman, 1997). However, little is known about the performance of these criteria in the case of the beta-mixture model. Our finding through both simulations and real examples is that neither criterion is suitable for the beta-mixture model. Both AIC and BIC seem to overestimate the number of components, leading to mixture models with excessive numbers of components. On the other hand, we find that the integrated classification likelihood BIC (ICL BIC) (Biernacki et al., 1998) is adequate, always selects the right model in our simulation studies and yields reasonable results when applied to real data.
-実行してみる
--単峰性入力をクラスタ数1,クラスタ数2で推定。
#ref(BetaMM1.png,,50%);
--ニ峰性入力をクラスタ数1,クラスタ数2で推定。
#ref(BetaMM2.png,,50%);
---icl.bicの値が小さいほうがより適切なモデルということ。クラスタ数の推定が正しくできている。分布の歪みにも対応できている。
-'''[[2012/08/14 ラボスピーチで紹介:http://shower.human.waseda.ac.jp/~m-kouki/pass/20120831_Riken_seminar_miyazawa.ppt]]'''
-'''[[2012/08/14 ラボスピーチで紹介:http://shower.human.waseda.ac.jp/~m-kouki/pass/20120831_Riken_seminar_miyazawa.ppt]]（リンク切れ）'''
#ref(BMM.png,,70%);
-2012/09/04 入力データそれぞれのクラスタ番号を求める・混合分布の各分布を別々にプロットする方法を、菊池研究室の浅井君がまとめてくれました。どうもありがとう。詳細は'''[[浅井君の報告:https://www.evernote.com/shard/s64/share/5041-s64-3f72f2e4b8b3bf901cc508edeed9bac3-1/?nbn=%E3%83%99%E3%83%BC%E3%82%BF%E6%B7%B7%E5%90%88%E5%88%86%E5%B8%83&dpn=asaitakuya&un=asaitakuya&source=sharedNotebookInvite]]'''を参照（限定公開です）
**Rの混合モデル関連リンク [#a9981854]
-[[CRAN Task View: Cluster Analysis & Finite Mixture Models:http://cran.r-project.org/web/views/Cluster.html]]
-[[苦労する遊び人の覚え書き/ガウス混合分布（ＧＭＭ）:http://qh73xe.jimdo.com/%E3%82%AF%E3%83%A9%E3%82%B9%E3%82%BF%E3%83%AA%E3%83%B3%E3%82%B0%E3%83%91%E3%83%83%E3%82%B1%E3%83%BC%E3%82%B8/%E3%82%AC%E3%82%A6%E3%82%B9%E6%B7%B7%E5%90%88%E5%88%86%E5%B8%83/]]
-'''[[ノンパラメトリックベイズ（無限混合分布）モデル:http://shower.human.waseda.ac.jp/~m-kouki/pukiwiki/index.php?%E3%83%8E%E3%83%B3%E3%83%91%E3%83%A9%E3%83%A1%E3%83%88%E3%83%AA%E3%83%83%E3%82%AF%E3%83%99%E3%82%A4%E3%82%BA%EF%BC%88%E7%84%A1%E9%99%90%E6%B7%B7%E5%90%88%E5%88%86%E5%B8%83%EF%BC%89%E3%83%A2%E3%83%87%E3%83%AB]]'''
-'''[[ノンパラメトリックベイズ（無限混合分布）モデル:http://speechresearch.fiw-web.net/restricted/index.php?cmd=read&page=%E3%83%A1%E3%83%AA%E3%83%BC%E3%83%A9%E3%83%B3%E3%83%89%E5%A4%A7%E5%85%B1%E5%90%8C%E7%A0%94%E7%A9%B6%2FInfinite%20Mixture%20Models&word=%E7%84%A1%E9%99%90%E6%B7%B7%E5%90%88%E5%88%86%E5%B8%83]]'''
