幾霜::残日録::2006/11/18 (土)

　昨日も書きましたが年代較正点を入れることができないようです。最終的に出力されるTime.treesで根までのNode Heightを1として枝長が出ているので適当に値を何倍かして計算しろということのようですね。しかしそれでは根点1点しか較正点にできないんですが・・・。そうか、PhyBayesが使われないのはそういうことなのか。使いにくくてもmultidivtime使うしかないのね。

追記 - 06:32:27
　ちゃんと計算すれば根点しか較正点にできないということはないな。しかし面倒だし、1点しか入れられないのは変わらない。これが勘違いだといいんですが。論文の方では複数の較正点入れているみたいですし。ドキュメント化されていないだけですかねぇ。

　1回1回の計算結果を見てみると、推定値の信頼区間はかなり広くなってますね。全くお話にならないくらいに。これはPL法と同様の傾向です。やっぱり大量の較正点を入れないと階層モデルに基づくMCMCはダメなんですかね。図はこれまでと異なり実線はPhyBayesによる推定値(事後分布の平均値)でエラーバーが95%信頼区間です。ベイズ推定した信頼区間とブートストラップによる信頼区間の違いはどういったものなのか、ちょっとまだよくわからないところがあります。真の値は昨日までのと同じです。

追記 - 06:22:27
　burn-inが少なかったので増やしたらどうかと思いましたが、大差ありませんでした。ありゃりゃ。

　とりあえず、シミュレーションを走らせ始めました。超疲れた。っつーかコントロールファイルのndataが使えればこんな苦労しないのに。ndataを指定してデータファイルに複数データセットを用意しておいても2つ目以降解析が行われない。おかげでディレクトリを分けたデータを1,000組用意する必要があった(ファイル名が指定不能な重要出力ファイルがあり上書きされてしまうため)。

　現在、burn-inは最初の200世代で、以降50世代毎にサンプリングで5,000世代計100系統樹をサンプリングするように設定して解析中です。しかもこれまで同様にデータ生成は異なる3つの分子進化速度変化パターンを仮定して行っているので上記の解析を計3,000回行うわけだ。すごい演算量だろうなぁ。あ、ちなみに速度の事前分布はexponentialです。これ以上はあまりやる気が起きませんねぇ・・・。

追記 - 04:40:29
　Linuxバイナリを使ってAthlon 64 X2マシンで2つずつ並列に解析を行っていますが、8OTU・2,000bpで約3分強かかってますね。3,000回の解析に9,000分ってことは150時間ですか。

追記 - 05:51:52
　burn-inが足りないかもしれないので猛烈にburn-inを増やした解析を一つ実行中。

追記 - 06:11:08
　burn-inが足りないっぽいので200から5,000へ増やして、3パターン各1,000回を100に削減。

追記 - 06:22:27
　burn-inを増やしても、あんまり結果変わらなかった・・・orz。

　期待通りの素晴らしい本です。理論分布図・度数分布図・散布図が豊富で分かりやすく図解されています。かつて苦労してイメージした画像がここにあります。

　ただ1点、PC用ソフトウェアとしてExcelを想定して書かれていることだけが残念でなりません。これがRを使っていたら最高だったんですが。Rの良いところはデータを生成するのが簡単なのと、実データが付属していて簡単にロードできるところなんで、その辺りは他のソフトでは今のところ太刀打ちできない。さらに統計教育用の作図デモ入りパッケージとかあったら最高なんですけどねぇ。とか言ってたら実はあったりして。