ぎょーむ日誌 2007-06-(01-10)

2007 年 06 月 01 日 (金)

• 0800 起床． コーヒー． 朝飯． 0925 自宅発． 晴． 0940 研究室着．
• 大雪山系稜線上の地衣類繁殖の階層ベイズモデル論文のつづき． まずは figure legends を完成させねば， と …… おうぷんおひすというか「わーぷろ」はつくづく使っていて楽しくない． いつも使っている Vim などテキストエディターであれば検索・置換・補完・画面分割などなど すごくすばやく簡単にできるのになぁ． あと， LaTeX のマクロその他の自動化・省力化機能あれこれは重要だと思うんだけど， わーぷろには実装できないんだろうね．
• とぶちぶち文句たれつつ， Results 節かきに着手する．
• いきなり調査シーズンになったためか， A 棟 8F の人口密度えらく低い． 1150 現在， 私ともうひとりいるだけ． まあ， 野外調査だけでなく， 今日は起学専攻の中間発表とかもあるんだけど． 聞くところによると この意味不明ぎみな名称を冠する 同専攻は常軌というものを必ずしも重視しない傾向があり， 院生全員の「中間発表の練習」までもが麗々しく挙行されただけでなく， 「練習」に専攻関係者全員出席が要求されたそーで …… そしてあとで聞いたところによると， どういう内容の発表であっても 「それは地球温暖化とどういう関係があるのかいちいち説明せよ!」 と命じられるそーで．
• Results 節もだいたい書けたので昼飯．
• 志水さんがファイルダウンロードできるディレクトリ準備し， そこに「同期」をとるための Makefile 書き．
• おっと References がまだだった． モデルまわりは引用文献数がむちゃくちゃに少ない． 太陽軌道三次元計算は高校数学・高校理科あたりの知識の 組みあわせで解決する問題なので特に引用すべき文献とかないんだけど， ベイズモデルまわりは本来ならもう少しいろいろと引くべきであった． しかしめんどうなので Clark and Gelfand (2006) だけですませてしまう．
• 何が「めんどう」かといえば BibTeX 使えぬこの「わーぷろ」 おうぷんおひす writer で手作業で References 作るってのは すごく苦痛で ……
• だいたいこの「わーぷろ」の紙面とゆーか画面とゆーかが真っ白 ってのはどうにかならんものかね …… そうか! なーんだ液晶ディスプレイの輝度 (明度?) おとせばよかったんだ! たしかに近ごろの LCD は昔とくらべて明るくなりましたねえ．
• とりあえずひととーりできたので印刷出力してひと休み． うーむ， やはり印刷出力しても醜い． とくに数式まわり．
• そして私の敵国語作文もヒドいものだなあ …… とひととーり修正してから， ディレクトリの同期をとって， 志水さんにメイル送信して …… 一件落着か? いやはやどうなることやら． LCD の輝度をもとにもどす． 時刻はすでに 1650．
• と脱力してから気づいたんだけど， (今回も，というべきか) MCMC 計算の収束判定については何も書かなかったな． library(coda) とかにいろいろな convergence diagnostic あったりするわけだが …… みょーな「検定」つかったものが多いわけで．

  検定では「違う」とは (確率論的に) 言えるけど 「同じである」とは言えない

  という基本をふまえると， 収束判定とやらには使いたくないわけで．
• ゐんばぐすは (複数チェインで計算してるときは) の推定計算結果を library(R2WinBUGS) で受け取ると， そこで Gelman-Rubin diagonistics やるんだよね． 下の図 (今回の地衣類繁殖の階層ベイズモデル; この図は今までぎょーむ日誌に置いてなかったことでもあるし) でいう R-hat (つまり ) が 1 に近ければ「収束」してます， と． おおざっぱにいって， chain 内分散の平均に比べて「全 chain まとめた」分散 (= chain 内分散 + chain 間分散) が何倍大きいか， とゆー数値．

• library(R2WinBUGS) の bugs() 関数が生成する bugs オブジェクト (上の図は plot(bugs オブジェクト) 図) lichen.bugs を print(lichen.bugs, digits.summary = 3) すると

  Inference for Bugs model at "/home/kubo/lichen/mcmc/lichen.bug", fit using
  WinBugs, 3 chains, each with 60000 iterations (first 30000 discarded),
  n.thin = 100 n.sims = 900 iterations saved
  
                mean     sd     2.5%      25%      50%      75%    97.5%  Rhat n.eff
  pbase[1]     1.446  0.237    1.013    1.263    1.448    1.608    1.925 1.041    66
  pbase[2]    -0.843  0.409   -1.635   -1.128   -0.840   -0.567   -0.041 1.043    62
  parea[1]     1.820  0.112    1.614    1.734    1.821    1.899    2.037 1.017   150
  parea[2]     1.535  0.162    1.226    1.424    1.532    1.647    1.856 1.034    76
  psri[1]      0.543  0.298   -0.043    0.342    0.539    0.752    1.087 1.034    77
  psri[2]      0.030  0.408   -0.749   -0.244    0.028    0.302    0.857 1.010   590
  iv.pface     1.291  0.288    0.822    1.104    1.247    1.450    1.930 1.012   240
  iv.pid       1.018  0.084    0.864    0.959    1.013    1.073    1.190 1.004   450
  pface[1]     0.018  0.408   -0.780   -0.245    0.012    0.305    0.808 1.029    71
  pface[2]     0.278  0.433   -0.506   -0.030    0.287    0.589    1.116 1.001   900
  ... (中略) ...
  pface[63]   -0.280  0.299   -0.812   -0.495   -0.303   -0.085    0.372 1.002   900
  deviance  3476.261 39.377 3400.000 3450.000 3473.500 3503.250 3556.000 1.000   900
  
  For each parameter, n.eff is a crude measure of effective sample size,
  and Rhat is the potential scale reduction factor (at convergence, Rhat=1).
  
  pD = 525.3 and DIC = 4001.5 (using the rule, pD = Dbar-Dhat)
  DIC is an estimate of expected predictive error (lower deviance is better).
  

  ここで Rhat についての簡単な説明もでるわけだが …… より詳しくは， たとえば R2WinBUGS 論文， つまり

  Sturtz,S., Ligges, U. and Gelman, A. 2005. R2WinBUGS: A package for running WinBUGS from R. Journal of Statistical Software 12:1-16.

  (Gelman web site のこの猛烈なる published papers からダウンロードできる; Gelman and Rubin (1992) の PDF file もある) なんかはどうかね， と眺めてみると，

  In this plot, the left column shows a quick summary of inference and convergence (R is close to 1.0 for all parameters, indicating good mixing of the three chains and thus approximate convergence); and the right column shows inferences for each set of parameters.

  としか説明してない． Gelman たちの Bayesian Data Analysis 本 (2003) の ``11.6 Inference and assessing convergence'' などに library(R2WinBUGS) のR-hat などに対応する定義が示されている． Gelman and Rubin (1992) の定義もこれとほとんど同じだけど， びみょーに異なるような気がする．
• また Gilks 本 Markov chain Monte Carlo in Practice (1996) の ``Chapter 8. Inference and monitoring convergence'' で Gelman が同じ内容を解説しているな． 先日の Clark ゾウ本にも ``7.5.1 Has It Converged?'' として Clark 御大の我流な導出が．
• 日本語の解説としては …… たとえば， 岩波・統計科学のフロンティア (12) 計算統計 II -- マルコフ連鎖モンテカルロ法とその周辺 -- の ``5.1 Gelman-Rubin の方法'' として種村さんによる解説がある．
• …… ということで， いまさらながらにして R-hat の正体が 少しは正確にわかったので， このあたりの記述も追加するか． R-hat とか持ち出さなくても昨日しめした MCMC trace 図で十分という気もするけど． トドマツ原稿とかにも追記せんといかんだろうな (こっちには MCMC trace 図がない)．
• などと書いてたらトドマツファイル一式が届く． ちらっとみて …… うーん， たいへんそう …… ということで本日は撤退． 2015 研究室発． 2030 帰宅． 体重 67.4kg． 晩飯．
• [今日の運動]
  • 腹筋運動 30 × 3 回． 腕立ふせ 10 × 3 回． スクワット 100 回．
  • ストレッチ
• [今日の食卓]
  • 朝 (0830): 米麦 0.6 合． 納豆． ワカメスープ．
  • 昼 (1320): 研究室お茶部屋． 食パン． チーズ．
  • 晩 (2150): 米麦 0.8 合． ネギ・ニンジン・ワカメ・高野豆腐の煮物．

2007 年 06 月 02 日 (土)

• 0700 起床． コーヒー． 朝飯．
• 洗濯． 怠業 …… としたいところだけど， ちょっと部屋を片づけないことには 生活にさしつかえるほどになってしまったので整理整頓 & 掃除． けっこうたいへんな仕事になった．

[この部屋にきて初めて]

ベッドの下のがらくたが (一時的ながら) 一掃された． 何か感銘のごときものを覚えてしまったので， 思わず写真撮影． ついでながらと言うか， 夜になってからこのベッドのグリル部を (がらくたの中から発掘された) 針金で補修． はい， 例の「愚行の柏」 事件で「不可逆的改造」しちゃったアレを 6 年後の今ごろになって ……

• 昼すぎまでかかって， あれこれ捨てるべきものの分別に成功． 当家にやたらとある， とわかったもの:
  • ズボンのたぐい: 先日の体重変遷図 (あるいはそれ以前の (なかなか) 15kg 減量に至らない記録) に示されているとーり， ここ 10 年間ほどの体重増減ハゲしく， それにともなって当然ながら 胴まわりの長さも伸縮繰り返してしまったがゆえに， さまざまなサイズのものが必要になったため， か?
  • 裁縫 (小) セット: がらくたにまみれてすぐに行方不明になるから
  などなど．
• しかし何とも汚い部屋に住んでたもんだな …… と， しばしぐったり脱力してから活動再開． 昼飯．
• Clark 御大 ゾウ本 の演習書・副読本とでもいうべき

  Statistical Computation for Environmental Sciences in R: Lab Manual for Models for Ecological Data

  が予想よりだいぶ早く届いた．

[ゾウ本 R & WinBUGS 篇]

表紙は「正」ゾウ本と同じ群ゾウ写真が使われている． 142 ペイジなのでウスい． amazon.co.jp では Models for Ecological Data (ペーパーバック) なる書名で売られている． 現時点で 2955 円 (私が買った時点では 2471 円だった)． 粕谷さん blog にもエントリーあり．

• あれこれながめてるうちに， 夜になってしまった． 晩飯．
• 当家は TV がないんでラジオきいてることが多いんだが …… 2130-2200 NHK ラジオ第二放送 NHK カルチャーアワー 文学の世界 の 4-6 月は 20 世紀イギリス小説 その豊かさを探る で毎回聴いてるわけでもないんだけど， 講師の小林章夫さんの語りがなかなか楽しめる． 何回か前の David Garnett の 狐になった人妻 は特にきあいが入ってました． 今回は George Orwell の 1984 年， 来週は Graham Greene の 第三の男 だそーです．
• 2320 自宅発北大構内走． 夜も晴れてるな． 2355 帰宅．
• [今日の運動]
  • 北大構内走 35 分間．
• [今日の食卓]
  • 朝 (0720): 米麦 0.5 合． ネギ・ニンジン・ワカメ・高野豆腐の煮物．
  • 昼 (1500): 米麦 0.6 合． 麻婆春雨．
  • 晩 (2015): 米麦 0.8 合． アスパラガス・ニンジン・ピーマン・ネギ・ショウガ・ シイタケ・イカの炒めもの． キュウリ・レタスのサラダ．

2007 年 06 月 03 日 (日)

• 0800 起床． コーヒー． 朝飯． 洗濯． 怠業 & またちょっと部屋の片づけしてから， 1020 自宅発． 今日もいい天気 (札幌の「五月」晴れ?)． 1035 研究室着． うだうだと web site まわりを少々ひねくる．

[イタヤカエデ巨木]

地環研前にある老齢 (?) の Acer mono (手前はよくわからぬマツ)． 駐車場の車との比較で大きさがわかると思います (樹高 15 m ぐらい?)． 今年も元気に展葉しているよーです．

• メイル書きあれこれ ……

R プログラミングに関して気づいたことなんですが，工学系・数学・物理学系
の分野ならいざしらず，おそらく他の多くの学問分野の R ユーザーたちにとっ
て，R が「第一プログラミング言語」 (最初に実用的に使うプログラミング言
語で，しかもこれ以外はほとんど使わない場合が多い) になるという状況が現
出しつつあるような気がします．

私が教えている範囲では，データ構造と制御構造を教えると，(とくに院生ぐら
いの若い人たちなら) わりとすんなりプログラミングできるようになるように
なるようです．ただし，プログラミングにおける関数という概念はなかなか使
いづらいようです (概念は理解できるが自分で独自の関数を定義できない)．複
雑な処理が必要になると copy & paste で似たような (しかし微妙に異なる)
コードを列挙していくことが多いですね．

いっぽうで，他の言語である程度プログラミング経験のある人は R ぐらい簡
単に自習できるはずだと取りくんでは見たものの，(やや独特なところのある)
R データ構造の操作で行き詰まったりしてるようですね．

• inline なる package があり， これは何と R コード中に C/C++ コードを書いとけば， それをコンパイル (R CMD SHLIB) して R 関数オブジェクトと関連づける， という作業を自動的にやってくれる …… という便利なのか邪悪なのかよくわからぬもの． Unix 系の OS の場合， たいてい R CMD SHLIB で問題なく外部モジュール作ってくれると思うんだけど， たいはんのゐんどーづ R ユーザーの環境ではダメでしょうなぁ ……
• また別の R ハナシで 間瀬さん のメイルのご指摘で気づいたんだけど …… R の配列 (というかたいていの場合は vector) の添字 (index) わざはだいたい「こうだろう」というふうに動くんだけど， 添字の範囲をはみだしたりしたときにちょっと変わった動作をする．

  > x <- 101:110 # 長さ 10 の integer vector
  > x[1] # 1 個めの要素
  [1] 101
  > x[1:3] # 1-3 個めの要素
  [1] 101 102 103
  > x[-1] # 1 個めの要素をのぞく
  [1] 102 103 104 105 106 107 108 109 110
  > x[-(1:3)] # 1-3 個めの要素をのぞく
  [1] 104 105 106 107 108 109 110
  
  > x[0] # じゃあ 0 個めの要素は?
  integer(0)
  > x[11] # 11 個めの要素は?
  [1] NA
  
  > x <- 99 # こんどは長さ 1 の vector 作る
  > x[1]
  [1] 99
  > x[-1] # x は integer かどうかわからないので
  numeric(0)
  > x[0] # これも同じ挙動 
  numeric(0)
  > x[2] # この挙動は長さ > 1 vector と同じ
  NA
  

  まあ x[0] としたときにともかく「長さゼロの vector (など)」 返してくるのでそれで問題ないわけだが．
• また別の R のハナシ． 私自身は特に必要性を感じてないんだけど …… R で作った図を「Illustrator みたいなそふと」 で修正・修飾してみたい， という需要はあるのかもしれない． で， いろいろな方法あるんだけど， 例によってどんな OS でも (つまり Linux でも MacOS X でもゐんどーづでも) 使える free software で解決する手段としては inkscape (Wikipedia 解説) を使うと良いのかもしれない． ただし R 側出力でも多少の工夫が必要であり， RjpWiki 内の Rのベクターデータを弄る人へ にあるように library(RSvgDevice) の devSVG() 使って

  	devSVG(file = "output.svg")
  	plot(...)
  	dev.off()
  

  というふうに SVG 出力して (SVG は XML な vector data file)， あとはこの出力ファイル output.svg を inkscape (下の図は Vine Linux 版) で open してやれば， 図をひねくりまわしたり文字あれこれを変えたりできる．

• 日曜日なのに， なンか院生密度が高いなあ (これは研究室にいながら推定可能) …… と， こわがっていて昼飯時間がおくれたわけじゃないけど昼飯．
• 地衣類繁殖の階層ベイズモデル， MCMC 収束に関してごく簡単に追記しておく． 志水さんはおいそがしいようなので， この論文化作業はひとまず suspend されるのかしらん?
• トドマツ原稿あれこれにとりくむ． まずは図表の修正． これはけっこう手間だ． そして本文なんだが …… うーむ．
• ということで質問メイルなど書いてから撤退． 1850 研究室発． 1905 帰宅． 晩飯の準備． 晩飯．
• [今日の運動]
  • うんどう休養日
• [今日の食卓]
  • 朝 (0830): 米麦 0.7 合． 麻婆春雨．
  • 昼 (1440): 研究室お茶部屋． 米麦 0.6 合． レトルトパウチドカレー．
  • 晩 (2100): 米麦 0.8 合． ホッケ． キャベツ・ニンジン・キュウリ・アスパラガスのサラダ．

2007 年 06 月 04 日 (月)

• 0800 起床． コーヒー． 朝飯． 0910 自宅発． 晴． 今日もいい天気． 0925 研究室着．
• トドマツ原稿あれこれ修正作業のつづき． まずは種子散布カーネル図まわりの修正． いままで， 「事後分布の (ホントに) てきとー値計算」 ですませていたところを
  • 79 本の母樹ごとの「代表散布カーネル」 (慎重に，というか) median() とったものに
  • 集団周辺散布カーネル (population marginal dispersal kernel) は「平均化」事前分布からの計算を (とりあえず) ヤメて， MCMC sample にもとづく積分 (つまり 79 母樹 × 300 MCMC step × triple chain 個ぶんの散布パラメーター使った平均値算出) にきりかえる
  …… で， とりあえず png() 出力してみて変更前と比較してみたんだけど， 「見た目」はそれほど変わってないようだ (左-右がそれぞれ旧-新， 新のほうは母樹ごとに曲線の色を gray() でびみょーに変えている)． 破線は marginal な種子散布， 「平均化」事前分布をつかおーが MCMC sample による積分だろーがほとんど同じように見える， ということで安堵．

• 図表まわりはこれでひととーりかたづいたはずなので 本文のほうの修正作業に．
• …… と取りくもうとしたんだけど， なぜか雑用メイルかきとか机のまわりの整理整頓に逃げてしまった．
• 逃げてばかりもいられないので， 作文なおし作業にもどる．
• あまりススまぬまま昼飯．
• 夕方までかかって作文修正と説明メイルかき …… ファイルアップロードしてメイル送信した時点で時刻はすでに 1720． うーむ， 時間かかるな ……
• また事務関係のわけわからん「ゑくせる」ファイル …… がんばって 30 分弱でオワらせる．
• そしてトドマツに関してはこれから 24 時間ほど suspend しそうな気配が …… うーむ， その間にイヤな査読が一個あるんだけど そっちを片づけてしまおうかな． 〆切は 10 日後なんだけど．
• で， 原稿よんでみたんだけど， やはりイヤになってしまったので撤退． 1815 研究室発． 1830 帰宅． 体重 67.4kg． 晩飯．
• [今日の運動]
  • 腹筋運動 30 × 3 回． 腕立ふせ 10 × 3 回． スクワット 100 回．
• [今日の食卓]
  • 朝 (0830): 食パン． キャベツ・ニンジン・キュウリのサラダ．
  • 昼 (1330): 研究室お茶部屋． 米麦 0.6 合． インスタントのヒジキ味噌汁． リンゴ．
  • 晩 (1910): 米麦 0.8 合． ニラ・ニンジン・ピーマン・シイタケ・干しエビ・ 卵の炒めもの．

2007 年 06 月 05 日 (火)

• 0830 起床． コーヒー． 朝飯． 1015 自宅発． 晴． 1030 研究室着．
• Occupancy Estimation and Modeling: Inferring Patterns and Dynamics of Species Occurrence …… という occurence つまり 「ある場所に A という species がいる・いない問題」に関する データ解析モデリング本がでてるな． amazon.co.jp で 詳しい目次 とか見られるようになっているんだけど， どうも空間相関とかきちんとは考えてない モデリングを紹介しているような気もする (内容みてないんで断定できないんだけど)． こういうたぐいのデータ解析って生態学会大会のポスター発表なんかでも よく見かけるんだけど， やってる皆さんは「これでぜんぜん OK」とか思ってるのかしらん?
• 値段たかめのような気がするんだけど， アリ科研費とかで買ってみるかな? アリの occurence とか調べないといけないようなハナシもあるわけだし．
• …… などと現実逃避ばかりもしてられないので， 昨日のつづき， すなわち読むと不幸な気分になる査読原稿にもどる．
• 論点をまとめるために原稿読みつつ査読報告書もぱらぱらと書きはじめてみる． あまりススまぬうちに昼飯．
• 昼飯後も投げ出したくなる衝動に耐へつつ査読報告かき …… 1658 送信． 一件落着． がっくり． 受け取ったほうもがっくりだろう．
• 査読途中の気分転換として， ポスターはがしその他なる雑用を処理． なぜか「北大祭」に地環研も巻き込まれてしまい， ポスター展示しなければならなくなったため．
• 地衣類繁殖モデルに関して， 志水さん担当部分の文章ファイルと Adobe Illustrator で描かれた大雪山地図が送られてきた． まずは地図をお茶部屋 iMac の Adobe CS の Illustrator で開く． 次にそれを SVG で保存して自分の ThinkPad (Linux 機) に転送してみる． SVG なら inkscape で開けるわけだが ……

• フォントまわりのびみょーなところに問題あるような． 要研究． 今日のところは Illustrator にもどって EPS で「別名保存」したファイルをつくり， それを Linux 機内 latex file 内で \includegraphics{} させる．
• 文章ファイルあれこれも OpenOffice & NeoOffice 方向に統合していく． お茶部屋 iMac に NeoOffice インストールしてみたんだけど， NeoOffice2.1.dmg が 136 MB もあってダウンロードもたいへん． まあ， さすがに OpenOffice ⇔ NeoOffice 間のファイルやりとりは何の問題もない …… いや， font がちょっとね．
• あれこれアップロードして志水さんにメイル連絡して， 1905 研究室発． 1920 帰宅 …… しかし忘れものその他が原因であちこち一時間ほど歩きまわることに． 2040 再帰宅 晩飯．
• 米がつきたので北 12 生協で買いもの． なぜか無洗米が安かった (10 kg 3180 円) ので買ってみる． 安い理由は， 混合比がきらら 397 : ななつぼし = 8 : 2 となる「純」道産無洗米， といったあたりか．
• [今日の運動]
  • うんどう休養日
• [今日の食卓]
  • 朝 (0930): 米麦 0.6 合． ネギ・シイタケ・豆腐・鶏肉のスープ．
  • 昼 (1320): 研究室お茶部屋． 食パン． リンゴ．
  • 晩 (2100): うどん．

2007 年 06 月 06 日 (水)

• 0740 起床． コーヒー． 朝飯． 0910 自宅発． 晴． 0925 研究室着．
• 昨日の夕方になって統合版をつくった地衣類繁殖の階層ベイズモデル論文の 原稿のみなおしのつづき． 主著者の志水さんとしてはこれを地衣類まにあのための Journal たる The Lichenologist に投稿したいらしい …… ざっと検索してみると， Bayesian な手法は地衣類まわりでは系統解析だけに使われていて 生態学的な問題ではまだ使われてないな (そういえば似たような検索を google scholar でもやったことあるな)． というか， そもそもいわゆる生態学的な問題をあつかった地衣類論文は そんなに多くないのである (分類・系統解析・地理分布の論文が多いみたいだ)．
• ということで地衣類モデル作文のなおし．
• 奇妙な表現のようだけど， Markov Chain Monte Carlo なチェインという意味で MCMC chains とゆー表現はけっこう使われてしまっている (私も使ってる)． まあ， 使われてしまっているからといって， これが適切かどうかよくわからんのだけど．
• Discussion 節もちょっと書いてみる． 書いてて 「そもそもこの二種類の地衣類が好む『明るさ』って違うんでは」 と今さら気になったので， 今さら図を作ってみる． まあ， このサンプル数では何とも言えないような …… あえて言えば Umbilicaria cylindrica (画像) は場所に好みがあるかもしれないけど， Lasallia pensylvanica (画像) はよくわからん． ついでにデータ構造も 岩(岩面(個体)) となってて ややこしーし． そう言や， 岩面の面積もわからないんだっけ．

• 作文なおし． いつまでたっても おうぷんおひす writer は楽しくない． テキストエディターの世界に帰りたいよー …… っていうか数式があろうがなかろうが LaTeX のほうが論文かきなんかはだんぜんラクだと思うんだが (頻出する学名のマクロ化とか)．
• MacOS X 使いの志水さんから (当方がアップロードしたファイルを) 「NeoOffice で open したら数式部分が□文字列挙状態」 ……?! いやはや． 当方 Linux 先方 MacOS X という OS 依存な問題かしらん?
• うーむ， Vine な OpenOffice の Math (oomath) の font を作業前に変えとく (default とする) ってのをやってなかったのがマズかったのかしらん? お茶部屋 iMac 上の NeoOffice では何の問題もなかったのになぁ …… と正しいのかどうかもよくわからぬ面倒な手作業を． こういう「おひす系そふと」なんぞ使ってしまうと， 人間にとってじつに屈辱的な労働を要求されてしまう状況が多発するね． まあ， 多くのユーザーにとってはそれが「安心感」なんだろうけどね …… でなきゃ， 皆さんかくも喜んでこんなソフトウェアつかうわけないもんなぁ．
• てな愚行に時間とられてしまってまた昼飯が遅れた． 昼飯．
• ばかげた単純労働もかたづいたので， 地衣類繁殖努力を推定する階層ベイズモデルの 事後分布をながめて議論すべきことを Discussion 節に書き込んでみる．
• ひととーり書いて， またファイルアップロード． 志水さんにメイル書き． 本日の地衣類まわりの作業はここまでか?
• 共著者のみなさんがおいそがしそうなので， トドマツ方面はしばらく動きがないもよう．
• ということで C 棟 7F に昨日はがしたポスター届ける往還 (この地環内の棟間移動は歩くとけっこうな高低差がある) にて， さて次は何をしようかと考える． ポスター配達そのものは， 受け取り担当の人がすでに帰宅していてからぶり．
• 私が何か関与していて目下のところ進行中の論文化作業は
  • 地衣類繁殖
  • トドマツ母子距離
  • シウリザクラ繁殖
  • アカマツは関係者の動きを静観中?
  といったところで (ホントはも一件あるけどコワすぎて書けない)， 次にやらねばならぬ作業の候補としては
  • アリ関連 の出口の見えないデータ解析 (deadline 8 月アタマ)
  • 屋久島葉寿命 の階層ベイズモデルの論文化
  あたりなんだろうが ……
• じつは他にも M2 問題というのがあって …… 趣味と実益と実業を兼ねて 修士論文のデータ解析に強制介入しているんだけど， 今年度は M2 は 5 人 (うん? もっといるのか?) もいて， しかも甲山さん M2 が 3 名も． 多種系における「重量分配」 推定問題はわりと単純な階層ベイズモデル化でなんとかなったわけだが (e.g. ぎょーむ日誌 2006-12-27, 2007-02-07) …… もっと他の形質 (つまり長さだの太さだのあれこれ) もあるんだよねえ． で， 憶測としてはこういった問題にはこれまでのような 一変量事前分布ではうまくいかない場合があって， 多変量事前分布， といってもおそらく多変量正規分布 + Wishart 分布ぐらいだろーけど， つかわんといかんのかもしれない， と．
• アリはかなりアブなくてとても嫌なかんぢだし， 屋久島葉寿命のハナシもさっさと終わらせたいんだけど …… M2 問題もやはり危険である． というのもパイプ樹木発表だの統計学授業だので 今年の 9-11 月はどたばたしそうな予感があって， ですね．
• 授業といえば昨日提出した事務書類ファイルの内容はこう． このまったく意味不明な内容はここに記録しといて， 腐れゑくせるファイルは廃棄．

• ともあれ， Clark 御大論文と著書でもながめて次なる計算にそなえて 勉強をすすめてみることにする． さすがに以前よんだときよりは理解できるなぁ ……
• 空腹になってきたので， 1925 研究室発． 1940 帰宅． 晩飯．
• [今日の運動]
  • うんどう休養日
• [今日の食卓]
  • 朝 (0830): 米麦 0.6 合． ネギ・シイタケ・豆腐の味噌汁．
  • 昼 (1440): 研究室お茶部屋． 食パン． リンゴ．
  • 晩 (2130): 米麦 0.8 合． ニラ・ショウガ・ニンニク・干しエビ・卵の炒飯． ネギ・シイタケ・豆腐の味噌汁． キャベツ・タマネギ・キュウリのサラダ．

2007 年 06 月 07 日 (木)

• 0800 起床． コーヒー． 朝飯． 0910 自宅発． 曇． 0920 研究室着．
• 0930-1040 地環研講堂で教授会． 会議といっても何かとくに議論があるわけでもなく， 書類の束が配布されてその内容を研究院長の岩熊さんその他が 読み上げていく， というもの． ということこで， 昨日のつづきで Clark 御大論文を読んでみた． いまながめてるのは， 1.5 年前の 論文紹介 でとりあげた Clark et al. (2003)．
• 以前みたときには二次元正規分布の事前分布 + Wishart 分布の超事前分布の 「つかわれかた」 というのがよくわからなかったんだけど …… 今日よみなおしてみてよーやく理解できたんだけど， この論文で紹介されてる階層ベイズモデルにおける 多変量正規分布の事前分布の分散-共分散行列は， biological (な解釈) には何の意味もなく， (極端ないいかたをすれば) ただ単に MCMC 計算の収束を速めるためだけに役だっているらしい， と理解できた．
• 複数の推定値のあいだにはどうしても (biological には何の意味もない) 「相関」ができてしまう (例: いわゆる直線回帰の「切片」「傾き」の推定値間の相関) ので， どうせなら事前分布にくみこんでしまえ (sampling の効率がよくなるだろう)， という発想だろう． ふーむ， やはり係数間の分散共分散行列を生成させても， こっちの「推定値間の相関」を見てしまうことになるよな ……
• つまりこういう状況だ． 個体 i で観測された何かの現象を説明する 統計モデル内に a + b x_i という線形部分 (linear predictor) があるとしよう． 説明変数 x_i は個体 i に固有なもので， 説明したい現象も一回だけ観測されたものとする (これは上の Clark 論文であつかってる例題と同じ状況)． このときにパラメーター a と b の両方に「個体差」をいれたモデリングをやりたいところなんだけど (たぶん現実には両方に「個体差」あるだろうから)， a + b x_i と足してしまうので両パラメーターに「個体差」いれてしまうと ややこしくなる (だけならよいが，しばしばわけわからなくなる)， と．
• で， Clark 御大方式 (と仮に呼ぶことにしよう; Bayesian modeling ではよくある手口) というのは a と b の事後分布がそれぞれ独立した事前分布から生成されるのではなく， パラメーター間の「相関」 (といっても解釈上の「biological な意味」とやらはとくにない) を考慮した二次元正規分布の事前分布から生成させる， と． 無情報な超事前分布から生成される (事後的な?) 事前分布は おそらく a と b のあいだに強い負の相関があるだろう (すると 2 パラメーターが実質的に 1 パラメーターに縮退する場合もある)．
• 私自身はこのような状況での階層ベイズモデリングにおいては， パラメーターと説明変数をよくみて， 「個体差」パラメーターの数を 可能なかぎり減少させる， という方法で対処してきた (これもよくある手口)． つまり上の例だと a には「個体差」あるとして b にはいれません， といったモデリングで面倒を回避できる． Clark et at. (2003) でもいくつかのパラメーターにはこの方策がとられている．
• Clark 方式と久保方式を比較してみるとこうなるだろう． 久保方式は多変量正規分布だの Wishart 分布だのもちださなくてすむぶん， まあ何となく「わかりやすい」という印象をもたれるかもしれない． しかしながらこの方式では状況がもっと複雑になった場合， たとえば観測が一回だけでなく複数回になったりすると うまくいかなかったりするだろう． これに対して Clark 方式はどういう場合でも柔軟に対応できそうなので， 拡張性のあるやりかたになっているのだろう．
• ところで， ここであつかってるパラメーター推定値間の相関に関して， biological (な解釈) には何の意味もない， と強調してるけど， べつにこれはそれが悪いとかいってるわけではない． 観測値間だの推定値間だのに生じる「相関」にはいろいろな種類があり， biological に意味のない「相関」もでてくるんだけど， 統計モデルのなかではそういうのも注意ぶかく扱わねばならないね， ということだ． 悪いのは「パラメーター間の相関」といえば 何でもかんでも 「とれーどおふだ!! 進化の帰結だ!! 多種共存だ!!」 とか無意味に騒ぎまくったりする理解の浅い一部の連中の所業である …… むろん狡猾なる Clark 御大はこんな簡単なところでしくじるはずもない． われわれ 職業的うそつき にとって用心ぶかさは美徳なのである．
• そして私をナヤませてるのは biological に意味が「ありそう」な 相関を多変量な事前分布で表現することなのだが， この論文にはそれに対するお手軽回答はない． 御大の近著論文 (まだ online early 版のみ) ``Resolving the biodiversity paradox'' はそのあたりちょっとふれているようなんだが ……
• ともあれ階層ベイズモデリングに興味あるヒトは Clark et al. (2003) を精読してみればよろしかろう， と思うわけです．
• Clark et al. (2003) ではいろいろと MCMC 計算わざが使われてるのだけど (Wishart 超事前分布のパラメーターの与えかた，など)， そのネタ本のひとつは

  Carlin P.C. & Louis T.A. 2000. Bayes and empirical Bayes methods for data analysis. Chapman & Hall / CRC

  …… で， この本もってるんだけど， どこにそういう Wishart 分布わざネタがあるのかよくわからん． 目次・さくいんではみつけられん．
• 見つからんはずで， もう一方のネタ本である Gilks et al. (1996) (Clark et al. (2003) では 1995 と誤記されてるが) のほうに掲載されてるじゃないか! p.307:

  ... To avoid this while still allowing the random effects a reasonable amount of freedom, we adopt the rule of thumb wherein ρ = n / 20 and R = ...

• 自宅冷蔵庫内の残りものかたづけるべく， 1240 研究室発． 走ると 7-8 分で帰宅できる． 昼飯． のんびりと皿洗いなどして， 自宅発． ちょっと雨． 1330 研究室もどる．
• トドマツ原稿に関して， 東大の練さんからいろいろとご指摘が． 修正作業にとりくむ．
• …… で思わず熱中してたら， 時刻はすでに 1415． 1400 から始まる別の会議には遅刻だな． ということで不参加決定． ちなみにこれは 「平成 19 年度第 1 回大学院環境科学院教授会」 とやらで， 午前中のは 「平成 19 年度第 3 回研究院教授会」 なのである． 組織を意味もなくひねくりまわすから会議が増殖するわけだ．
• しかしこれって「定足数」とかいう条件くりあーしないと， 「再試合」になったりするわけだな． 会議中は集中して勉強できることだし， あまりサボらないようにしよう．
• 1450 ごろトドマツ原稿の修正作業終了． やっぱテキストエディター (vim) & LaTeX & R 作図を組みあわせた環境は修正がラクでいいわ …… ファイルアップロード & 練さん・後藤さんにメイル送信．
• 論文原稿といえば， 志水さんからの連絡で， 当方の腐れおうぷんおひす数式エディタでフォント設定を改めると， 志水さん Mac の NeoOffice でも問題なく数式が表示されるとわかった．
• 教訓: Vine Linux が packaging した openoffice.org の数式エディタの font 設定は 「他の OS でも見れそうな font」 が default となるようあらかじめ設定変更しておく (すくなくとも sazanami のたぐいはダメみたい)
• 昨日はからぶりに終わった 「地環研内はげしく上下移動ポスター提出の旅」 は本日は問題なく終了． ついでに気分転換の運動になる．
• 当節の学術ぢゃーなるにしては発行頻度が低い Environmental and Ecological Statistics (← まあこういうマニアな方向性なので 論文あまりないのかもしれないけど)， 最新号の目次がきた． ふーむ， ``Simulating correlated count data'' (DOI URL: http://dx.doi.org/10.1007/s10651-007-0008-1) といったおもしろそうな内容のものあるな． 相関のある Poisson 分布由来のカウントデータだしたければ， まあふつーはこの平均値の対数を多変量正規分布にしてやる， といったワザがよく使われるんだけど， この方式では
  • 必ず overidspersion があるものと仮定せねばならぬ
  • 平均値が低い場合に低い相関しか許容されない
  といった欠点あるので， それを改善する方法を考えました， というもの． まにあだ． たしかに私は重要な問題だと認識するんだけど ……
• げ， やはり恐れていたとーり， いきなり 蟻類研究会大会 の要旨だせとかいういつものごとき唐突なる大統領命令が発せられたよ …… 来週月曜日提出， か． なーんかでっちあげないといけないね．
• とゆーことで， また岩倉さんのエゾアカヤマアリ修論などひっぱりだしてみる． うーむ ……
• 関係ないけど， ゐきぺでぃあに アルゼンチンアリ の項目が新しくできたな．
• 1835 研究室発． 1850 帰宅． ちょっと時間おそいけど洗濯機まわしてみる． 1930 洗濯終了． 晩飯．
• 社会保険庁 の 年金個人情報提供サービス とやらで， アカウント申請するといま話題の個人年金データをながめることができる (ただし共済年金に関しては支払回数わかるだけで詳細はわからない)． 私の場合， 国民年金 → 厚生年金 (いまは亡き NASDA 傭兵の独房群時代) → 共済年金， と変化してきたわけで …… ふーむ， 国民年金に関しては「納付済」で 110ヵ月がうまっているな．

• [今日の運動]
  • またうんどう休養日 ……
• [今日の食卓]
  • 朝 (0840): 米麦 0.6 合． キャベツ・タマネギ・キュウリのサラダ． シシトウ炒めもの． ブナシメジの味噌汁．
  • 昼 (1255): 米麦 0.6 合． ニラ雑炊． キャベツ・タマネギ・キュウリのサラダ． シシトウ炒めもの．
  • 晩 (2030): 米麦 0.8 合． タマネギ・卵の炒めもの． キャベツ・タマネギ・キュウリのサラダ．

2007 年 06 月 08 日 (金)

• 0700 起床． ねむい． コーヒー． 朝飯． 0845 自宅発． 晴． 0900 研究室着．
• また石狩浜のエゾアカヤマアリデータながめる．
• そして志水さんと大雪山系の地衣類繁殖の階層ベイズモデル 論文かきの面談． また図などをいろいろと追加することに． 午前中はこれでほぼ終了．
• ECOAS 方面のメイル書きしてから昼飯．
• トドマツ修正原稿があちこちから届いたので， 昼飯後はそのあたりをまとめる作業． ファイル upload & メイル送信， で 1400 ごろひとまず終了． その後にまたちょっと Table まわりを修正．
• …… さて， 次だけど …… 各個撃破を企図して， 地衣類論文の追加作図にとりくんでみるか．
• しかしなぜかしらその企図を皆さんに察知されてしまって (いったいどうやって!?) 各方面からのメイル十字砲火の撃ちすくめられて身うごきままならぬ状況に． だれか助けてー
• なンかメイル送信がどの SMTP server 使ってもおかしいような …… 北大の腐れ mail gateway とやらがまた腐れているのかしらん?
• ひとまず返信かいて …… しかし難問である ECOAS 方面の八甲田の湿原のハナシ (お， 松山大会サイト内に 講演要旨 が) はすぐに返信かけないので， 論文を熟読することに． 熟読しても難しいものは難しい．
• そしてトドマツ方面はどうなることやら …… と皆さんのどたばたの推移を傍観者的に見物してみる． うーむ， だいじょうぶなのかしらん． そして， とりあえずトドマツ方面はしばらく動きがないだろう， と楽観してみる．
• そしてふたたび八甲田湿原群の 割算値の絶望の闇 の中にずぶずぶと沈んでいく …… 数十分ほど沈んでいくと， 闇の中にヒトすじの光が見えてきた． これは [効率] = [光獲得量] / [重量] といった割算値問題ではなく， [光獲得量] = [重量] × [係数] なる問題に他ならない (ホントに重量に比例するのか疑わしいところだけど， この [効率] とやらの趣旨はそもそもそういうことだ)． 例によって例のごとく， [重量] はあたかも offset 項的な部分ね．
• さらに， だ． これは群落内で観測された (とゆーか憶測計算された) 資源つかみどり問題なので， [光獲得量] は植物種ごとに独立ではなく， 自分がぶんどれば他人のとりぶんは減り， 邪魔なお隣のさばればあちきの獲得量削りとられてしまう， という， まあディリクレ分布的な状況ね …… 形式的にはむしろ Log normal- Multinomial logistic 的な状況とでもいいますか?
• ようするに， こちらの研究室の院生たちと研究をすすめてる樹木個体内の 重量分配モデル と同じたぐいの統計モデリングだ． 組みこまねばならぬ (fixed and random) 要因おおくて イヤになるけど …… イヤかどうかはともかくとして， 計算に工夫が必要だな．
• 問題は割算値がひとつだけなら上のごとくでよろしいわけだが …… じつはすごくまぢめに考えると， これって [重量] が原因で [光獲得量] が結果であると決めつけてしまう 光ぶんどりあいのハナシでもないのかもしれないね …… 世の中， まぢめに考えてはいけない問題だらけね．
• 1945 研究室発． 2000 帰宅． 体重 67.0kg． 八甲田割算値やせ， だな． 晩飯の準備． 晩飯．
• [今日の運動]
  • 腹筋運動 30 × 3 回． 腕立ふせ 10 × 3 回． スクワット 100 回．
• [今日の食卓]
  • 朝 (0800): 米麦 0.7 合． ネギ・ブナシメジ・卵の炒めもの．
  • 昼 (1220): 研究室お茶部屋． 食パン． リンゴ．
  • 晩 (2050): 米麦 0.8 合． ショウガ・ニンイク・ホタテのバター炒め． コマツナのゴマあえ．

2007 年 06 月 09 日 (土)

• 1000 起床． コーヒー． 朝飯． 洗濯． 怠業．
• 今日もいろいろとメイルが …… こういう査読依頼， どうしたものかしらん?
• また八甲田湿原群の絶望割算値とらっぷから脱出する問題の検討， とか． 複雑なデータ構造をもつ問題ほど， 割算値をつかった場合にとりかえしがつかないほど解析がめちゃくちゃになる， という気がしてきた． つまり， 複雑な状況においてこそ 素朴な基本わざの組みあわせで解決する方策を模索せよ， ということか …… で， 当代においては Bayes なやりくちは 基本わざたちの「接着剤」となり， コンビネイションわざを発動させることできます， と．
• 1815 自宅発． 1830 研究室着．
• 八甲田湿原群の脱割算値メイルかき． 作文に苦慮する．
• すっかり遅くなってしまった． 2305 研究室発． 2320 帰宅． 晩飯．
• [今日の運動]
  • うんどう休養日
• [今日の食卓]
  • 朝 (1100): イソップベイカリーの巨大クロワッサン．
  • 昼 (1500): イソップベイカリーのバゲット． レトルトパウチドのスープ (「緑の野菜と岩塩スープ」がおいしかった) キャベツ・タマネギ・キュウリのサラダ． プラムトマト．
  • 晩 (2350): イソップベイカリーの雑穀パン． タマネギ・ニンジン・ネギ・豆腐のスープ． キャベツ・タマネギ・キュウリのサラダ． プラムトマト． コマツナのゴマあえ．

2007 年 06 月 10 日 (日)

• 0840 起床． コーヒー． 朝飯． 洗濯． 怠業．
• JavaScript コード書きで遊んでしまった …… たとえば， テキストボックスにカーソルをもっていく こちらで紹介されてる ようなやりかたでできる， とわかった． なかなか奇妙な言語よのう ……
• 1150 自宅発． 晴． 札幌駅周辺うろうろ． 寝るとき用に 「T シャツを長袖にしたようなもの」 あったら便利かも， とさがしてみたんだけど， この季節そんなものは売ってるわけない， と判明した． また本屋とかにも寄り道してしまって， 1340 研究室着． だんだん雲ってきた．
• メイル書きしてから昼飯．
• どーしたもんかなー， と思いつつも …… 査読依頼謝絶メイルかき． あまりこういうことはやらないんだけど， やむをえない理由で．
• 大統領命令で アリ研究会大会 の要旨を明日提出せねばならぬことに …… この要旨とやらをでっちあげるにしても， このときに「みょーに自信ありげな作文」 をものするためには解析するデータを多少なりともさわっておくのが 有効だろうと思われたので …… 岩倉さんの石狩浜エゾアカヤマアリのデータファイルの解毒 (脱ゑくせる化) に苦闘してみる． とりあえず tab 区切りテキストとして export． 「セル内コンマ区切り」 などという解毒難易度向上わざが使われてるためだ． もちろんこのままでは R とかでは読めないので， Perl で変換フィルターを作っていく． ケガれ言語たる Perl を使う理由は， セル内にうめこまれた不規則な日本語文字コメントをうまく処理せねばならぬ， という事情があるため ……
• 外はいつのまにか雨．
• どういう形式のデータ構造に変換しようかいろいろとナヤんでいて， ふと気づいたんだが …… うーむ， 多項 logistic モデルにおいて 「同時に複数の応答」 ってのは …… なんだか許されるような気がしてきた． つまりサイコロでいえば 「1 と 3 が同時にでる」 といった状況 …… しかし， よくわからんけど， 「1 から 6 までが同時にでる」は許容されない? ともあれ， これってけっこうスゴいような． まあ， (サイコロの例でわかるように) ふつーそんなカタチで現象を記録しないんだけどね．
• …… そうそう， このデータの中で， 拉致られアリは「一回だけ」逆襲に転じていたんだよな …… この野外実験はいろいろなところからアリを拉致してきて， それを「敵地」コロニー近傍にほうりこんで， 侵入者としてあつかわれるて どういうメにあうか (地元ワーカーと接触 5 回分) 観察しましょう， というデータ． 拉致られアリは不安なのか基本的に「逃げ腰」で， 噛みつかれることはあっても相手を攻撃することはない …… しかし， 330 個体 × 5 回の観察のなかで， 一回だけ「逆上」してるんだよね．
• 残念ながらこの逆襲現象はモデル化のアイデアがないので， ここに記録だけして変換後のファイルからはとりあえず削除しておこう．
• とりあえず R で読めるようにしてみた．

> head(d)
  antid   colony     home    hometype tr30m trial aggdisplay allogrooming
1     1 hoshioki hoshioki supercolony     0     1          0            0
2     1 hoshioki hoshioki supercolony     0     2          0            0
3     1 hoshioki hoshioki supercolony     0     3          0            0
4     1 hoshioki hoshioki supercolony     0     4          0            0
5     1 hoshioki hoshioki supercolony     0     5          0            0
6     2 hoshioki hoshioki supercolony     0     1          0            0
  antennation biting ignorance others trophallaxis
1           1      0         0      0            0
2           0      0         1      0            0
3           0      0         1      0            0
4           0      0         1      0            0
5           0      0         1      0            0
6           0      0         1      0            0

> summary(as.factor(apply(d[,7:ncol(d)], 1, sum)))
   1    2    3 
1611   38    1 

• 全 1650 回中， 「一回の観察中に複数の行動が観察された」 は 39 回だけ， か．
• ばててきたので撤退． 1930 研究室発． 1945 帰宅． 晩飯の買いもの． 晩飯の準備． 晩飯．
• [今日の運動]
  • 腹筋運動 30 × 3 回． 腕立ふせ 10 × 3 回． スクワット 100 回．
• [今日の食卓]
  • 朝 (0930): イソップのバゲット． イソップ パン生活つづく． タマネギ・ニンジン・ネギ・豆腐のスープ． プラムトマト．
  • 昼 (1410): 研究室お茶部屋． イソップのバゲット． タマネギ・ニンジン・ネギ・豆腐のスープ．
  • 晩 (2100): 米麦 0.8 合． ニラ・ピーマン・ニンニク・卵炒飯． キャベツ・ニンジン・キュウリのサラダ． プラムトマト．