KuboLog 2007-05-23

苦情・お叱りは，たいへんお手数かけて恐縮ですが，久保 (kubo@ees.hokudai.ac.jp) までお知らせください．

2007 年 05 月 23 日 (水)

0750 起床．コーヒー．朝飯． 0930 自宅発．晴． 0940 研究室着．

さいきんの Ecology Letters は毎号のように Bayesian & MCMC な論文が掲載されてるな …… これもおもしろそうだ．

Subhash R. Lele, Brian Dennis, Frithjof Lutscher. 2007. Data cloning: easy maximum likelihood estimation for complex ecological models using Bayesian Markov chain Monte Carlo methods. Ecology Letters (OnlineEarly Articles) 10:???-???. doi:10.1111/j.1461-0248.2007.01047.x

数ヶ月前に，とある (やや工学より) proceedings の査読やったら Amazon ギフト券 3000 円もらってしまった．さすがは「やや工学より」，カネもちだねえ．ということで，さっそく Bayesian data analysis 本購入にぽーんとつかってしまった．

トドマツ原稿にとりくむ．そうだな，まあ discussion 節の改定文案かきから始めるか．

えーと，当方愛用のテキストエディター Vim で dos な (Shift-JIS) なファイルきめうちで開くには …… vim "+e ++enc=cp932 ファイル名"，と．

いつものごとく遅々たる作文を 1230 ごろまで．さーて，作文はいったんおいて，また事後分布まわりの計算をやらんといかんね．まあ，これがひかえめにいっても面倒きわまりないもので．

いきなり 30 分ほど 7F に拉致されて，甲山さん・山中さんの Global COE ヒアリング練習につきあわされることに．明日が本番だそうで．

ちょっと計算すすめてから昼飯．

また散布距離分布問題． R で (事後分布から MCMC sampling された標本の) ``quantile'' 値だす方法いろいろありそうだけど，おそらく一番簡単には， quantile() 使う方法で，たとえば 95% 区間はこう計算される．

> x <- rnorm(100000, 0, 1) # ここでは例題としててきとーな乱数
> quantile(x, probs = c(0.025, 0.975))
   2.5%   97.5% 
-1.9662  1.9530

さて，「平均散布距離」とやらの 95% 区間だけど …… 79 母樹それぞれに 600 個の (事後分布からの) 標本値があるので，「事後」散布カーネルは合計 47400 個あるわけだ．とうぜん，それぞれが平均値 (sqrt(0.5 * pi) * sigma) もつわけで，その中央値 (50% 値) と 95% 区間はこうなる．

normal2D(sd0)           33.6 (28.2, 41.3)
normal2D(sd0 + sdMt)    31.3 (7.2, 210.6)

上の normal2D(sd0) が (LMM ではなく GLMM 的な意味あいでの) ``random effects'' ぬきモデルで (sdMT <- 0; 蛇足ながらこいつの平均値は 34.0 m)， normal2D(sd0 + sdMt) は 79 母樹差 random effects ありモデル (蛇足ながらこいつの平均値は 51.6 m)．後者は，いやはや，ナイスな混合モデル的ばらつきかたですね．

とびこみ R こんさる．「Y 軸のラベルを水平にしたい」ですか．しかもよりによって「ぼうぐらふ」とは …… ! ともあれ回答例送信．

x <- 1:10 names(x) <- LETTERS[x] par(las = 1) barplot(x, horiz = TRUE)

[とりあえず help(par()) を]

par() の help はなかなかごちゃごちゃしてるんだけど， horizontal とかで検索すればすぐに par(las = ...) にいきあたります，と．

トドマツ原稿 discussion 節かきなおしつづく．

現実逃避というわけでもないけど， (院生密度がなぜかひくいので) お茶部屋で亀山さんややのんびりとツガザクラ雑談．いろいろ発見が続いてますねえ …… ところで，雑種 F1 はできるけど F2 ができない，という現象を説明するには複数遺伝子座の相互作用を考えないといけないのだろうか?

またトドマツの paper/section/discussion.tex の修正作業などなどを続けて …… まあ，本日はこんなところかと make pdf && make save の rsync でバックアップして，「こう変えてみました」メイルかいて送信．時刻はすでに 1855．

1940 研究室発． 1955 帰宅．体重 67.6kg．

晩飯後に原稿にとりくんでみたんだが …… 仕事すすまん ……

[今日の運動]

うんどう休養日

[今日の食卓]

朝 (0840): 米麦 0.6 合．納豆．ナメコの味噌汁．
昼 (1400): 研究室お茶部屋．食パン．チーズ．リンゴ．
晩 (2030): 米麦 0.8 合．ニラ・エノキダケ・卵の炒めもの．

ぎょーむ日誌 2007-05-23

2007 年 05 月 23 日 (水)