ぎょーむ日誌 2004-12-(21-31)

2004 年 12 月 21 日 (火)

• 0810 起床． 朝飯． 昨晩からのつづきで空間相関問題の検討に没頭する …… かなり整理できてきた． コーヒー． 0950 自宅発． 曇． 1000 研究室着．
• とりあえず， 面妖ぷろじぇくと plotnet の 意味不明予算消化陰謀のかたぼうをかつぐ発注あれこれを 北大生協 に． さて．
• お茶部屋では朝から加賀百万石の大吟醸酒の tasting が ……
• 昨日から考えている空間統計学問題， じつは一生懸命かんがえたあげくに， pseudo likelihood method を再発明しただけでは， という気がしてきた． いやいや EM アルゴリズムの再発明か? ともあれ， すでに誰かが似たようなことを考えているのは， 間違いなさそうな気がする．
• と勉強してるといきなり部屋の電源おとされる． 院生部屋の局所的停電 (電力消費過大が原因) なおそうとして配電盤いじってたら， こっちまで落してくれたそうで …… 試行錯誤して停電解除． 教訓．
  • 8F 配電盤のカギは A701 などといったタグがついている
  • 問題の部屋のブレイカーはいったん 完全に off にしてから再度 on
• 昼飯．
• 1500 より 講座セミナー ， 本日は塩寺さんで， まだあまり高くないイタヤカエデの葉面積などは， 何で決まるか， というもの． いつもながら， シュートレヴェルの解析は難解である． 特に， シュート内の葉の相互作用． このあたり， きちんと解析するのは難しいだろう． ということで， 問題はシュート内の葉の相違をどこまで表現するか， に帰着される．
• セミナー後も葉雑談． さらに「高くないイタヤカエデ」といっても， 樹高 30-300cm ぐらいの範囲にあるので， 明るさ環境が同等かどうかよくわからない． というのも， 300cm になるまでに生きのびた個体ならば， 良い環境にいる可能性が無視できないからだ．
• 1815 研究室発． すすきのの居酒屋 月夜のたぬき で忘年会． まあ， さわがしいというかなんというか …… 一次会終了後に撤退． 2150 帰宅．
• [今日の運動]
  • うんどう休養日．
• [今日の食卓]
  • 朝 (0830): 米麦 0.8 合． ニラ・モヤシ・ショウガ・鶏レバの炒めもの．
  • 昼 (1340): 研究室お茶部屋． 米麦 0.8 合． 切り干しダイコン・コンニャク・サバブシの煮物．
  • 晩 (1910): 忘年会．

2004 年 12 月 22 日 (水)

• 0830 起床． 朝飯． コーヒー． なんとなくシュート内資源分割問題の検討に没頭してしまう． どうもまっとうな計算は無理だろうな． 1000 自宅発． 晴． 1010 研究室着．
• 昨日の停電事件の余波． 一年ほどまえに 甲山さんが衝動買いした PowerMac G5 は 180W の無停電装置では 停電時に支えきれない， と判明した． というのも CPU たる PowerPC G5 一個で 50W 以上の電力を食っているからだ． PowerMac の場合， 巨大な液晶ディスプレイの電力供給も本体から行なわれるし． かとーさんが登校してきたので， かとー iMac で実験してみると …… これは 180W でなんとか持久できるようだ． ということで， こちらにつけかえる．
• 甲山さんが研究室見学の学生さんと昼飯に行くので， 私もなんとなくその一行に加わる． ふーむ， 百周年記念会館 1F の「きゃら亭」がやけに混んでると思ったら， 本日はランチバイキングだったのか．
• 昼飯後にその学生さんあいてに講座概要説明し， あとの詳しい説明は院生の皆さんにおまかせしてしまう．
• また plotnet 予算消化雑用． あやしげ GIS ソフトウェア (私は使わないだろうけど) 購入で目標額達成をもくろむわけだが， いろいろ書類かかされたり， ですね． これだから独占ソフトウェアというやつは．
• ネット途絶 …… 原因をさぐってみると， A 棟 3F の建物 hub が原因らしい， ということで強制リセット． 原始的だ． これが有効だったかどうかわからんけど， ともかくその直後に復旧した．
• A-802 院生部屋電力途絶． また修繕． PowerMac G5 が電力くいすぎてんのかしらん， と思ったけど院生各自で使ってる電気ヒーター (出力 800W ぐらい) に問題あるようで． もっと抜本的な問題は， この A 棟 8F のつくりが極寒地北海道にあるまじき， やすっぽいできで隙間風が入ってくることで …… 甲山さんと相談する．
• R プログラミングこんさるてぃんぐ． 今日のお客さんは淡々と修論を片づけつつある平林さん． gene duplication のごとき copy & paste プログラムだったので， ループを使えばと書き直してみる． しかしあとで考えてみると， プログラミングに不慣れなヒトたちにとって， ループのような制御構造と， 処理の関数化のどちらのほうが理解しやすかったのだろうか?
• また研究とは関係ない脇道にそれてしまった …… swftools といういろいろなフォーマットの画像を Flash に変換するもの を見つけたんで， インストールして遊んでしまった (configure したら avifile がない， と警告されたのでそれもインストールしてみた)． まあ， 現時点ではぱらぱらまんがの生成ぐらいしか， 使い道を思いつかんけど．

• 夕方からぎょーむススまず， 統計学勉強ばかりやってしまった． またしても …… 2120 研究室発． 2130 帰宅． 運動． 晩飯．
• ここ数日， 日付変更時刻すぎても ぎぶす分布問題が脳内 working メモリから消せない． アタマの中で空間相関のある模様がマルコフ連鎖で明滅 …… この問題を具体的に調査できるようになったら， そのときに検討すべき手順としては
  1. Metropolis-Hastings アルゴリズムで autologistic な Gibbs 分布を生成してみる
  2. Hubbell や Jorge たちの (自称?) ``Gibbs sampler'' で 1. と同じ性質 (クリーク密度などなど) をもつ分布が生成されるか確認する
  3. 2. のやりかたで Gibbs 分布にしたがうパターン を生成しつつ推定値を改善していく --- iterative ``最尤推定'' (らしきもの) をやってみる
  4. maximum pseudo likelihood (MPL) 法で得られた推定値の生成する パターンと比較する
  5. これまでの MPL 推定値の ``補正'' 法の結果とも比較してみる
  …… といったあたりかなぁ． はたして 3. でやっている推定計算は， MPL 法のありがちな variation に過ぎないのか， EM アルゴリズムの不出来な亜流なのか， それともやってみたらわけのわからない推定値しか得られぬ 単なるごみなのか ……
• [今日の運動]
  • エアロバイク 45 分間．
• [今日の食卓]
  • 朝 (0840): 米麦 0.5 合． ニラ・モヤシ・ショウガ・鶏レバの炒めもの．
  • 昼 (1220): きゃら亭． またしても 甲山さんにごちそうになってしまいました．
  • 晩 (2300): 米麦 1.0 合． 切り干しダイコン・コンニャク・サバブシの煮物． ハクサイ・エノキダケ・煮干の味噌汁．

2004 年 12 月 23 日 (木)

• 1030 起床． うーむ． 朝飯． コーヒー．
• またまた空間集中性の推定問題にとらっぷされる …… かなり極端な観測データが得られたときに， いま検討してる推定法で破綻なく推定できるかどうか考えてみる． とくに問題なくいくような気がする． Gibbs sampler が「うねり」続けてくれるならば．
• 1600 自宅発． 雪． 考えながらふらふらと歩き回る． 1650 研究室着．
• 苫小牧 M2 平尾君との空間集中性メイルやりとり返信を書いていると， アタマの中が整理されてきた． いつもながら， 大学院生には教えられるねえ …… 平尾君の整理によると， 空間相関モデリングには 単なる autoregressive model と conditional なやつがある， ということになる． 前者の autoregressive model ってのは私の理解の中では maximum pseudo likelihood (MPL) 法に分類されるものだ (言われてみると， たしかにどの文献の数式表記でも MPL になってるな …… ああいうのって仮の表記ではなく， ホントにそのまま推定計算やるつもりだったのか?)． いっぽうで， 言うところの conditional autoregressive model はベイズ推定とのこと …… ここがポイントだよねぇ．
• 私の考えてる計算方法はベイズ推定というより， むしろ混合モデルの MCMC 推定版というかんぢだ． 言うまでもなく， 空間パターンそのものが Gibbs 分布にしたがう確率変数である， という扱いになるわけで …… 平尾君メイルの整理でアタマがすっきりしたのは， べいぢあんなやりかたで確率分布化されるのは 推定されるべきパラメーターだけでなく (私はこっちだけだろうと思ってた)， 必ずしも推定対象とは言えない空間パターン集合も含む， といったあたりが， じつにようやくにして理解できた， というところだ．
• これらいっさいがっさいを逐次的に推定していくのは， EM アルゴリズム的である， と言ってしまってよいのかもしれない．
• そうかそうか …… 上では「必ずしも推定対象と言えない」 と書いた， Gibbs 分布にしたがう空間パターンそのものを推定してみせるのが， Bayes 推定 + MCMC 法あわせわざによる 「画像修復」 問題というやつになるわけだ． つまり空間パターン集合のうち MCMC 法で得られた中で 「最良」 のものを選びだすのが (というか「探索する」というべきか) 最大事後確率 (MAP) 推定 …… うーむ， いろいろわかってきたぞ． いまさらながら．
• 確率分布を表現する個々の変数・パラメーターのうち， 直接観測できた数量以外は， なんでもかんでも確率分布から得られた確率変数とおいてしまう Baysian 流やり口は知識だけでなく実感としてわかった． 同時に， もしこの Baysian 推定計算の考えかたが受け入れられているものであるならば， いわばそれを簡単化している混合モデル MCMC 法も通用しないはずはあるまい． まじめな検討の対象になりうるような気がする．
• ベイズなモデルにしろ， Gibbs 分布から得られたパターンを「参照」する 混合モデルにしろ， 注意すべきは Gibbs sampler の設計と， その意味づけだな． 少なくとも， 生態学における 空間集中林冠ギャップ形成問題 と 空間集中な leafminer 寄生ハチ問題については， Gibbs sampler は簡単なものになりそうだし， ランダム変量を生成する Gibbs 分布の意味も明確であるように見える． このあたりはここまでしつこく考え続けてきた問題だ．
• というふうに， 11 月はじめの鹿児島研究集会いらい勉強にかまけてきた成果あれこれが， やはりよーやくにしてというべきか， アタマの中で相互の関係を明らかにしつつ使いものになるようになってきた， というかんぢだ．
• さらにネット上で文献あさる統計学勉強がつづいてしまう …… ほらほらやっぱり， 経済学の分野とかでは Bayes であっても空間相関の部分は pseudo likelihood 的なまま， というような論文もあるぞ． 何でもいいから 推定値らしきモノえられればそれでヨシってのも， なんだかねえ ……
• 統計学的手法ユーザーとして， pseudo likelihood method の最もまずいと思えるところは， 推定したモデルをもちいて乱数を生成できないところだ (もちろん間瀬さんの林冠動態論文で比較対象用に計算されていたように， MCMC に MPL 推定値をほうりこめばパターンは生成されうるけれど …… これは推定モデルと生成モデルが一致していない)． 「乱数」 が生成できない確率論的モデルは不完全というか， 危険すぎてパラメーター推定には使う気になれない． 同じ理由で quasi likelihood なモデルもイヤなんだよね ……
• 今日も事実上の怠業状態ではあった． どうなることやら． 2140 研究室発． あいかわらず， 雪． 2150 帰宅． 運動． 晩飯．
• 空間相関問題について， 一定の解決というか整理整頓ができたので， 今夜はアタマの中で Gibbs 過程から生成される パターンが明滅しなくなった …… しかし生活周期ずれたまま． さっさと寝なくては．
• [今日の運動]
  • エアロバイク 55 分間．
• [今日の食卓]
  • 朝 (1050): 米麦 0.7 合． ハクサイ・エノキダケ・煮干の味噌汁．
  • 昼 (1530): 米麦 0.7 合． ハクサイ・エノキダケ・煮干の味噌汁．
  • 晩 (2350): 米麦 1.0 合． ハクサイ・エノキダケ・煮干の味噌汁．

2004 年 12 月 24 日 (金)

• 0825 起床． 朝飯． コーヒー． 今朝は統計学問題にそれほど長時間はとらっぷされずに 0930 自宅発． 晴． 0940 研究室着．
• 来週月曜日に今年最後 (おそらく) の 苫小牧 もうでに行くことになりました．
• 三浦さんにその書類をおわたしするついでに， 甲山さん思いついてしまった plotnet 怪しげ消化案に関して， 交渉をお願いする． ちょい難しいかもと思いつつ，
  「……というような口実で， なんとかならないでしょうか?」
  「…… 私，イヤです ……」
  「ええっー?」 (硬直)
  「いえ，久保せんせーが，ではなくて． 会計掛はなんでも『だめ』としか言わないんで」
  この予算消化作業のいんちきさと不安感を端的に表現する会話である．
• 空間相関問題から「釈放」されたので， ようやくにしてぱいぷ樹木原稿なおしにとりくむ．
• ひととーり指摘点を原稿内でチェックしていって昼飯．
• 昼飯後もぱいぷ樹木原稿．
• ちょっといきづまったので吹雪の北大構内を散歩．
• なぜか A802 電力問題をしらべる電機屋が本日やってくる． 配電盤には予備の回線があるようなので， これを使って A802 のコンセントを増設， 負荷分散を試みてみることに． ただし工事は年があけてから．
• 会計掛問題， 1 月中旬の 監査を突破できるんなら勝手に買いやがれ (知ったことぢゃねーよ)， というような投げやり返事をもらったという． たぶんこの問題に関しては監査うけても問題ないだろうという気がする． しかしながら， なんというか， こういう連中が研究の流れを制動しうる， というのは当の会計掛も含めて 誰にとっても不幸なことだというべきだろう．
• 院生のノート PC がネットにつながらなくなる． このテの事故はひさしぶりだな， と思いつついくつか試験してみる …… 本体くみこみの LAN カードまわりがおかしい， と特定． またかよ， Fujitsu (そしてカニマークの ethernet チップ) …… さいわいにして， PC カード式の LAN カードが私の机のひきだしの中にあり， その Fujitsu 機には PC カードスロットあったんで， 使ってみる． XP 上で WNT 用ドライヴァーが動作したんで (もはや数年前に製造中止になったブツなんで XP 用ドライヴァーが入手できなかった)， XP の阿呆なコントロールパネルと格闘して バイパス設置に成功．
• 修正作業はあまり進捗せづ 1950 研究室発． 雪． 2000 帰宅． 体重 73.8kg． 運動． 晩飯．
• 局所明るさ問題， すこし面白い解決策を思いつく． 感度解析をやればよいのだ． 一石二鳥．
• [今日の運動]
  • エアロバイク 50 分間．
• [今日の食卓]
  • 朝 (0950): うどん．
  • 昼 (1330): 研究室お茶部屋． 米麦 0.8 合． 飯しかないので， お茶部屋に食品棚に残留していた 粉末中華スープで湯づけ．
  • 晩 (2200): 米麦 1.0 合． チンゲンサイ・ピーマン・油揚・卵の炒めもの．

2004 年 12 月 25 日 (土)

• 1000 起床． 朝飯． コーヒー． 怠業． 昼飯． 1550 自宅発． 曇． 1600 研究室着．
• またぱいぷ原稿にとりくむ．
• 院生部屋 A802 の防寒性をたかめる作業． ブラインドのついてない窓があったんで， お茶部屋から移植してみる． さて， これで少しはマシになるのか?
• で， ぎょーむのほうの作業は進捗せづ …… すごすごと撤退． 2120 研究室発． 2130 帰宅． 体重 74.2kg． 運動． 晩飯．
• 現実逃避的隠者的思考として， 「なンか， しち面倒な集団遺伝学問題は， だいたいにおいて MCMC 法つかった (よく考えると別に MCMC でなくてもいいのだが) 『隠れマルコフモデル』 で説明できてしまうのではないかしらん …… 世の中なんでも隠されている!」 といったアイデアに脳内メモリーを占拠されてしまう． メモリー管理が意のままにならずメモリーリークしたままになったりするのは， 人間とダメ OS の共通点である． いやいや， ここしばらく大学院生たちの隠匿度のたかいセミナーばかり 拝聴してきたせいか?
• 蛇足ではあるが， 集団遺伝学の難しい問題に対して MCMC が常に有効なわけではない． というのも， 観測データの標本もとである隠れ局所集団はちょっと油断すると 「尤度ゼロ」 状態 (統計力学用語でいえばエネルギー無限大) に落ちてしまうからだ． むろん「ええ， ここで， ここで奇跡のよーな突然変異が発生したのです」 と言いぬけるのは常道であるが， クリスマスの夜も終わりつつあることだし， おそらく世の中それほど奇跡にあふれてはいまい (そもそも冬はそれほど遺伝子流動ないだろうし)． 他にも隠れ集団サイズがどーのこーの， といった難点もある．
• しかし， 連鎖している多遺伝子座のランダムサンプルを生成してみよ， といった問題には MCMC はひどく便利だな …… と， いまさらながらに気づく．
• [今日の運動]
  • エアロバイク 55 分間．
  • 腹筋運動 30 × 3 回． 腕立ふせ 5 × 3 回．
• [今日の食卓]
  • 朝 (1010): 米麦 0.5 合． わかめスープ．
  • 昼 (1500): ロールパン．
  • 晩 (2330): 米麦 1.0 合． コマツナ・油揚煮物．

2004 年 12 月 26 日 (日)

• 0930 起床． 朝飯． コーヒー． 怠業． いやはや．
• 昼飯くってから 1520 自宅発． 晴． 札幌ではとっくに真冬日 (最高気温が氷点下) の日々に突入している． 1530 研究室着．
• うーむ， またいつものごとく， Gibbs 分布 / MCMC 法 / 空間相関 / 生態学・遺伝学への応用あれこれ 問題にトラップされてしまった．
• Gibbs 分布をランダム変量とみる混合モデルってのは， 少しヘンてこかもという懸念はあるけれど …… 形式的・意味論的には問題なさそうな気がする． ある空間パターンの出現尤度で加重すべきか? …… Gibbs sampler つかうんだから， そんなことは考えなくてもいいような気がするんだけど．
• 思考暴走とまらず， ついつい (必要もないのに) 連続空間問題との対応なんぞも考えてしまう …… こりゃいかん， 現時点での私の理解では， まったく対応がつきそうにない． 考えるのヤメ．
• アタマわやくちゃになりそうな多遺伝子座混合モデルなんだけど …… とにかく基本は 「全遺伝子座ペアの出現確率」 を計算し， そこからの逸脱でもって連鎖不平衡を定義すりゃあいいってこと． ペアを押さえれば， 三重以上の連鎖もそこに含んでいるわけだから． 連鎖という現象である以上， ねじれた相関はありえないはずだ．
• 晩飯． 2130 研究室発． 2140 帰宅． ハンダごてと太い針金でエアロバイクの壊れペダルを修繕する． 安物はいかんねえ． 運動．
• バスだと朝の札幌市内突破は時間かかりそうなんで， 明日は JR で苫小牧に行くか． 0810 発だと 0912 苫小牧着， と．
• [今日の運動]
  • エアロバイク 55 分間．
• [今日の食卓]
  • 朝 (1030): 蕎麦． コマツナ・油揚煮物．
  • 昼 (1450): うどん． コマツナ・油揚煮物．
  • 晩 (2040): 研究室お茶部屋． 平林さんにめぐんでもらったスパゲッティー． レトルトパウチドのカレー．

2004 年 12 月 27 日 (月)

• 0700 起床． 朝飯． コーヒー． 0750 自宅発． 大雪． 0800 JR 札幌駅着． 0810 快速エアポートに乗って同発． 0826 北広島で苫小牧行き普通に乗り換え． 0912 苫小牧着． ここまで 1410 円． 今回はタクシーに乗って 0930 苫小牧研究林 (TOEF) 着．
• 正午すぎまで， M2 平尾君と村上さんと leafminer 寄生蜂モデリング相談． お， すでに glmmML による解析まで． さてさて， 寄生率そのものの空間相関はどうすれば表現できるだろうか? モデリングの前段階， としてだ …… アイデアにつまったんで， 寄生個体間の距離の平均 vs 分散の関係と， 寄生個体の場所を randomized してみたパターンを比較してみれば， というような． はたして， これは使いものになるのか?
• 院生の皆さんと村上さんにお供して昼飯を食べに出る． ひさしぶりに藤そばの蕎麦を食べることできた．
• 午後から夜まで鍋嶋さんと苫小牧直径モデリング相談 …… 今回は短く終るかも， と思っていたんだけど， 私の説明がだらだらしていたり， とちゅーから健太さんの混合モデル質問介入あったり …… ともあれ， 私にとってもアタマの中の歳末大整理になりました．
  • TeX2Word つかわねばならんかもしれん …… (LaTeX 使いこなせそうなヒトに LaTeX おススめできぬ現状がつらい)
  • deviance 計算あれこれ必要かも …… めんどくさいから full model の尤度を 1 にしてしまうか?
  • sensitivity より elasticity? …… しかし elasticity って固有値の解析， という特定の状況でもっとも便利になるような気もする
  • ランダム変量の確率分布 (gamma distribution) 描いてみる必要あるな …… 樹木の個性が指数分布的になってたりして
• また平尾君・村上さんと議論． 平尾君はすでにくだんの平均 vs 分散 plot 作っている． なんたる早わざ． で， 私の反省なんだけど， これでたしかにランダムパターンからずれてるのはよくわかるんだけど， 「集中斑」 の個数によって 「平均小 & 分散小」だの「平均大 & 分散大」 だののパターンがあって， これだけ見ていてもよくわからんかも． 「平均小 & 分散小」だったら集中斑はひとつ， といった憶測はできるけど．
• うーむ， 樹木ごとの寄生率なる割算値の「地図」 を作らんといかん， ということかしらん?
• 最後は M2 平松さんの glm() in R こんさるてぃんぐ …… 苫小牧でも R がいよいよ普及しつつあり， めでたいことである． さて， 今回の問題は glm() が計算を「なげて」しまうこと．
• 調べてみてわかったこと: モデルが複雑すぎて glm() の行列計算がとまっていた． えー， glm() 愛好する大学院生のみなさん， いくらデータがたくさんあっても (そして指導教官の無茶な「やれ」指示があったとしても) いきなり glm(y ~ x1 * x2 * x3 * x4 * x5 * x6, ...) というような 「ちょー多項 (多重?) 交互作用」 いりモデルを計算させてはいけません． まずは glm(y ~ x1 + x2 + x3 + x4 + x5 + x6, ...) あたりでモデル選択などやってみて， 「現象の説明に使えそう」 な項だけを (もしホントにホントに必要だというなら) 交互作用にしてください． 交互作用とやらはせいぜい二項ぐらいまで．
• 平尾君の車で送っていただく． 2030 TOEF 発． 2044 JR 苫小牧駅着． 2048 普通に乗って同発． 2119 千歳で快速エアポートに乗り換え． 2151 札幌着． 2210 帰宅． 晩飯．
• [今日の運動]
  • うんどう休養日．
• [今日の食卓]
  • 朝 (0720): 米麦 0.7 合． コマツナ・油揚煮物．
  • 昼 (1240): 苫小牧市音羽町 「藤そば」 大ざる 880 円．
  • 晩 (2240): 米麦 1.0 合． チンゲンサイ・モヤシ・ニラ・卵の炒めもの．

2004 年 12 月 28 日 (火)

• 0820 起床． 朝飯． コーヒー． 0950 自宅発． 晴． 1000 研究室着．
• お， 生態学会の 地区会 の案内がいつのまにか．
• 昨日， 苫小牧の平尾君におしえてもらった 「R 用のまたべつの混合モデル計算関数」 を調べてみる． これは Jim Lindsey (R 開発にもいろいろ貢献した有名人らしい) が開発したけれど， なぜか CRAN には入っていない package repeated (Nonlinear Regression and Repeated Measurements Libraries) である (Lindsey の R code ペイジ からダウンロードできる)． repeated だけでなくそこらにある一連の package をまとめてインストールしたほうがよい． example が相互参照したりしてるんで．
• うーむ， glmm だけでなく gnlmm --- 'gnlmm' fits user-specified nonlinear regression equations to one or both parameters of the common one and two parameter distributions ... --- なんてモノまで含まれているのか．
• glmm から出力される結果オブジェクトをみると， sd なる coefficient が含まれている． これは overdispersion というかランダム変量のばらつきをあらわす パラメーターである． また， glmm の結果オブジェクトは stepAIC() への入力として使えるようだ．
• これって glmmML より便利そうだな …… しかしもっと慎重に調査してみなくては， と調査にかまけてしまった．
• 同じデータを与えて推定計算させてみると， 数値計算法がちがうせいか glmmML とはびみょーに結果が異るな (私が調べた数値例については， AIC はだいたい一致していた)． 現時点で判明している glmmML() より glmm が優れている点:
  • ロジスティック回帰で glmm(cbind(r, n - r) ~ x, ...) わざが使える (というか family = binomial の場合にはこれしか受けつけんようだ)
  • (先ほども書いたけど) stepAIC(glmm(...)) できる; summary(glmm(...)) もできる
• おお， repeated には隠れマルコフモデル (hidden Markov model; HMM) の推定計算関数 hidden() まであるのか． 「隠れ」好きなヒトはどーぞ．
• などというようなことばかりにかまけてしまった． 昼飯．
• いやー， なんというのか …… 「お勉強もーど」 からなかなか抜け出せない状態が続いてしまう． こういう擬似無限るーぷなのです:
  1. 統計学の教科書や論文をよむ
  2. すべては理解できないけれど， 書かれていることの一部は自分たちの問題解決に かなり有用そうだとわかる
  3. 実際に応用するにはどうしたらよいか， と定式化を工夫する
  4. すると文献中よくわからなかったあたりも 徐々に理解できるようになり， その確認のため さらに別の文献なんぞも読んでみて 理解度を試験してみたくなる
  5. 2. にもどる
  これにはまりこむと， なんとなく「勉強してる・仕事してる」感があって， 他のことが手につかなくなるんだよね． 研究商売やってる以上はこういう「お勉強」は必要不可欠なんだけど， 世間とは無縁ではありえないぎょーむ上の観点からすると， 「悪循環」としか言いようがない．
• 問題は 3.→4. の過程で …… 「えーい，お勉強なんぞにかまけてる場合ではない」 と文献類はどこかにしまいこんでしまっても， アタマの中では 3. の検討が止まらないのである． 地図をずっと見つづけていなくても山歩きは継続できる， というような状態だ． で， アタマ内検討がうまくいってもいかなくても， 数日以内には「地図」上での現在位置確認の欲求やみがたく， また「るーぷ」がまわりはじめるのである．
• てなことで午後も終わってしまい …… 1920 研究室発． 1930 帰宅． 運動． 晩飯．
• [今日の運動]
  • エアロバイク 45 分間．
• [今日の食卓]
  • 朝 (0900): シリアル．
  • 昼 (1330): 研究室お茶部屋． 米麦 0.7 合． チンゲンサイ・モヤシ・ニラ・卵の炒めもの．
  • 晩 (2230): 米麦 1.0 合． チンゲンサイ・ニラ・エノキダケの炒めもの．

2004 年 12 月 29 日 (水)

• 0940 起床． 昨晩は 2330 すぎに寝たのに， 起きるのがおそい． お勉強とやらで壊れたアタマの中の修復に時間がかかるのか? 朝飯． コーヒー．
• 純然たる怠業， ではさすがにいごこちが悪いので， 洗濯しながら怠業． あまり変わらんか．
• 1600 自宅発． 曇． 1630 研究室着．
• さすがに時間がなくなってきたんで， ぱいぷ樹木原稿なおしに着手．
• 甲山さんは年末になると， 年賀状をどたばたと印刷するためだけに登校してくる． それは良いとして， なぜか出産直後の雪野さんもそれにまきこまれている (嗚呼)． そして， すでに出産を終えたのに， 雪野さんお得意の (そして私を困惑させる) 妊婦ジョークをとばす． 年賀状あて名書きソフトウェアは， MacOS X の腐れた Classic 環境で動くので とらぶるが多発する． したがって， かとーさんもまきこまれる． 甲山さんと私はそれを傍観する．
• 甲山さんとぱいぷ樹木原稿に関する相談すこしばかり．
  • P_max 1.5 倍増問題のめんどうは 甲山さんにおしつける
  • NPP 倍増問題は呼吸量ほぼ一定説で説明する (これは甲山さん発案)
  • 光量増減問題については 「一次元問題ほどお手軽な世界ぢゃないんですよ」 と開きなおってみせる
  ふーむ， 相談すると展開が速い …… 場合もある．
• まんぜんと検索してるうちに思い出した． そういやむかし発生学で「動物」の体軸について 習ったんだよな …… 動物では
  • 前後軸 (anterior-posterior axis) = 頭尾軸 (cepahlic-caudal axis)
  • 背腹軸 (dorsal-ventral axis)
  • 左右軸 (left-right axis)
  • 先端基部軸 (proximal-distal axis)
  となる． ところで植物では 「枝の根元-先端方向の『軸』」 を proximal-distal とよんでいいのだろうか? あまりこういう言いかた見かけないよな．
• ちなみに私は root-terminal とよんでいた． これもあまり見かけぬ言いかただ．
• うだうだと作業してたのでまだまだ終了には遠い …… 2330 研究室発． 2350 帰宅． 晩飯．
• [今日の運動]
  • 腹筋運動 30 × 3 回． 腕立ふせ 5 × 3 回．
• [今日の食卓]
  • 朝 (0950): シリアル．
  • 昼 (1430): 米麦 0.9 合． チンゲンサイ・ニラ・エノキダケの炒めもの．
  • 晩 (2400): うどん． レトルトパウチドカレー．

2004 年 12 月 30 日 (木)

• 0830 起床． 朝飯． コーヒー．
• 苫小牧 の平尾君からのメイルで教えてもらった， geoRglm package で午前中がツブれてしまった． こりゃー， おもしろいおもしろい (平尾君は寄生率の再推定結果を送ってくれたんだけど， そちらは放置してしまう)．
• 概要は doc/geoRglmintro.pdf (Unix 系なら /usr/lib/R/library/geoRglm/doc/geoRglmintro.pdf) とか作者 Ole Christensen とか， あるいはお手軽には R News の短い解説文:

  Christensen, O. F. and Ribeiro Jr, P. J. (2002). geoRglm - a package for generalised linear spatial models. R News, 2(2), 26-28. ISSN 1609-3631.

  が良いだろう．
• geoRglm で使ってる generalized linear spatial model (GLSM) に関する R News 解説文の核心部はわずか十数行に要約されるわけだが， このところアタマの中がすっかり 空間相関 / MCMC / 面倒推定問題 / べいづ な状態になってる私は (めづらしいことに) 即独即解してしまった …… ひとことで乱暴に説明するなら， つまりは各地点で定義される線形予測子 (linear predictor) が「すっごーー … (果てしなく長い) … ーい 多次元正規分布」 になる一般化線形混合モデル， ってコトでございますよ． スっとんでいるけどうまい解決策だなぁ． MCMC との相性も絶妙だし …… というのも多次元正規分布ならば 条件つき確率分布の計算は 単なる行列計算に帰着できるからだ (分散・共分散行列の逆行列を計算したりするだけ) …… って自力で思いつかなかった私がアホなのかもしれないけど．
• この「線形予測子が (あたかも) ランダム変量」 なる久保流説明はヘンタイ的という他ないけれど， 形式的にはこれで間違っていないし， 自分としては一番おちつきが良い． まあ， もうちっと普通にいうと， 「各地点に空間相関をもつ隠れランダム変量 (ランダム変量は常に隠れてるわけだが)， これらを線形予測子に含む 一般化線形混合モデルをべいぢあん推定する」 というあたりか．
• ああ， べいぢあん， か． やっぱり そっちのほうが推定とかラクかもね． こういうのは．
• そしていま気づいたんだが， 上の乱暴説明にあるよーな線形予測子全体を 多次元正規分布に放りこむのは愚策というべきで， ランダム変量だけ MCMC するほうがよい． そうすればその平均はゼロとおけるからだ (分散・共分散行列だけ考えればよい)．
• 私がここしばらく あれこれ考えてる方式とはびみょーに異るところが興味ぶかい．

(項目)	久保式	GLSM
ランダム変量	ぎぶす分布にしたがうランダム空間パターン (離散値) そのもの	多次元正規乱数の要素 (連続値) を各地点にばらまく (まあこっちもギブス分布にしたがうわけだが)
推定方式	Bayes 推定もできるけれど， どちらかというと EM アルゴリズムつかった 一般化線形混合モデルの 最尤推定をねらっている	(おそらく) Bayes 推定のみ
尤度は?	あるパラメーターセットのもとで決まる 空間パターン集合と観測パターンのあいだの 確率論的距離の期待値	べいぢあんな事後分布にしたがう確率変数
Gibbs Sampler	なんとなく全地点一斉更新方式だけを考えていた (各地点を逐次的更新方式も使えるけど)	各地点逐次更新方式だけを考えている (たぶん …… マニュアルとか詳しく読まんとわからんけど; 全地点一斉更新も原理的には可能)
MCMC めんどくさ度	かなりラク	一回だけとはいえ， 莫大な個数のちょー巨大逆行列計算が必要??
説明 めんどくさ度	少し口ごもる …… かも (しかし GLSM 方式の説明をそのまま借用すれば)	かなりラクなはず …… 面倒をすべて「隠れ変数」 に押しつけることができるんで

• じゃっどん， GLSM 「ちょー巨大逆行列」問題， これは空間相関の及ぶ範囲を限定してしまわば， いづくんぞ難点あらんや， いやない (反語) …… しかもほぼすべての地点で「使いまわせて」おトク感ありますわよ! …… 境界条件うまく処理すれば全地点で使いまわせる．
• で， べいぢあんな推定やったら MAP 推定値とか計算せんといかんのだろうな …… アタマでは「最尤推定値と同じだ」 とわかってはいるんだけど， なンかみょーな違和感あってですね．
• 1130 自宅発． 晴． 1140 研究室着． いかん， 怠業状態．
• 昼飯． 昼飯後に上の GLSM うだうだを書きだしたら， またまた数時間をついやしてしまった．
• う． そしてかとーオフィス内完全怠業状態 3 時間． 年末ぼけだ．
• ちかん研 のわけわからん改組 「環境科学院」 にふりまわされる学部生から 誤爆的な問い合わせが …… とりあえず 新しい専攻科構造の表 (実質上ほとんど意味ナシ) をよく見て自分が連絡すべき相手をさがしてくださいと連絡する．
• で， 今朝がた送られてきた平尾君の R プログラムを拝見する …… う． すっげー copy & paste coding． 標本数 3 ながら 「いまどきの M2 は copy & paste coding を好む」 と結論してよいのかもしれない． まー， データ解析用プログラムを自分で書けるって， それだけでたいへん立派だけどね (たとえば 私の同世代なんかはそういうことできる生態学研究者はほぼゼロ)．
• ともあれ送っていただいた報告メイルみると， Jim Lindsey の repeated package の glmm() 関数は推定の安定性において glmmML() よりマシとは言えない， といった貴重な内容も． ふーむ …… 自作関数で Nelder-Mead 法とかで問題を解かせるか， 小標本のときは混合モデル使わない， というような対策が必要だねえ．
• なンか密度依存性とかは log() でいれたほうがよいのかもしれんねえ． その前に計算プログラムをどうにかしたいんだが …… というような返信．
• 霜月さん と BIC 談義． 私の調べたこと:
  • ふつーの BIC の定義: BIC = -2 × (最大化対数尤度) + (パラメーター数) × log(標本数)
  • Schwarz (1978) のモデル選択規準: SC = (最大化対数尤度) - 0.5 × (パラメーター数) × log(標本数)

  Schwarz. 1978. Estimating the dimension of a model. Annals of Statistics, 6: 461-464.

  • BIC は最小となるモデルを選び， SC は最大となるモデルを選ぶ; BIC = -2 SC
  私の憶測:
  1. (この項は事実) Schwarz (1978) で AIC との比較を論じつつ SC を提案する (このときには BIC という名前はない)
  2. だれかが 「AIC と比較しやすいように」 というような理由で BIC = -2 SC を定義した
  3. BIC と SC は「実体は同じ」という説明が， いつのまにか「形式的にも同じ」 という誤解が一部にひろまった
  4. まちがいが文献にのる (1) 「BIC = SC」
  5. まちがいが文献にのる (2) 「BIC = (最大化対数尤度) - 0.5 × (パラメーター数) × log(標本数)」
  (後記: まちがい (2) としている 「BIC = (最大化対数尤度) - 0.5 × (パラメーター数) × log(標本数)」 という定義は複数の教科書でみつかる …… もしかしたらまちがいでないのかも; 翌日のぎょーむ日誌参照)
• 2000 研究室発． 2010 帰宅． 運動． 晩飯．
• [今日の運動]
  • エアロバイク 45 分間．
  • 腹筋運動 30 × 3 回． 腕立ふせ 5 × 3 回．
• [今日の食卓]
  • 朝 (0840): シリアル．
  • 昼 (1230): 研究室お茶部屋． 「北欧」エビカツパン．
  • 晩 (2230): 米麦 1.0 合． キャベツ・ブナシメジ・コンニャク・油揚・煮干の味噌汁．

2004 年 12 月 31 日 (金)

• 0730 起床． 朝飯． コーヒー． 0920 自宅発． 晴． 0930 研究室着．
• BIC 「ふつー定義」調査のつづき． もし， ものごとが多数決で決まるなら

  BIC = -2 × (最大化対数尤度) + (パラメーター数) × log(標本数)

  が「ふつー」である， と言っても良さそうだ (e.g. R の lme4 や nlme に含まれる BIC() の定義，あるいは多くの教科書)． ところで Ripley 御大の著書 ``Pattern recoginition and neural networks'' の glossary には

  BIC has two similar meanings. Akaike (1977, 1978) introduced `information criterion B'. Schwarz (1978) introduced something which has become konwn as a `Bayesian information criterion'. Although most references mean Schwarz's BIC, to avoid confusion this is also know as SBC (`Schwarz Bayes Criterion'). Both penalize the deviance by log n times the number of p of free parameters for n examples, but Akaike's has O(p) terms not depending on n.

  などとある． Schwarz's BIC は deviance と p log n なる本質さえおさえていれば形式は何でもヨシなのかもしれぬ． ``Akaike (1977, 1978)'' は入手困難な文献で， 何が書いてあったか不明．
• BIC のマイナー (?) な定義， つまり Schwarz (1978) の (まだ BIC と呼ばれてなかった) 規準 そのまま

  BIC = (最大化対数尤度) - 0.5 × (パラメーター数) × log(標本数)

  を掲載している教科書の例として， ``Learning Bayesian Networks'' (Neapolitan, 2004)， あるいは ``分子系統学'' (長谷川 & 岸野, 1996) など (後者は 霜月 さんにご教示いただいた)．
• とまあ， こういう現状なので， 「ここで用いる BIC の定義は○○であり， これを [最大 | 最小] 化するモデルを選択する」 と明記しておけば何でも OK なんかねえ …… AIC と BIC を「比較」してみたければ ふつー定義のほうがぜんぜん便利だし， 比較なんぞ意味ナシで Schwarz (1978) の最初の定義だけが唯一無二のものだ と考えるならまいなー定義つかうしかない．
• 蛇足ながら Rissanen (1989) の MDL (minimum description length; 最小記述長) 基準 なるものがあり (「基準」 (standard? index?) と「規準」 (criterion?) の使いわけのきじゅんがわからん)，

  MDL = - (最大化対数尤度) + 0.5 × (パラメーター数) × log(標本数)

  これは Schwarz's BIC と同等である．
• もひとつ蛇足で， ややこしいことに Akaike (1980) 提唱による ABIC (Akaike's Bayesian Information Criterion) なるものまである． これは階層ベイズモデル (ここでは経験ベイズモデルだけを考慮してる?) の最大化対数周辺尤度と超パラメーター (hyper parameter) の個数できまる規準である …… うーむ， 苫小牧直径成長モデルのモデル選択規準は現在のままで OK， ということだと理解してよいのだろうか?
• お， なぜか A 棟のスチーム暖房が稼働している．
• なんつーか， いやいやながらというようなまいんどせっとで， ぱいぷ樹木原稿のなおしにとりくみはじめる． 難しくはないけれど， めんどうだ．
• 昼飯． そのへんに誰かいるよーないないよーな．
• 夜までかかって， ぱいぷ樹木「水まわり」の補足説明かき． つまりは， C++ で長々と書かれてるのを (てきとーに省略しつつ) 敵国語にすればいいだけのことで．
• 1920 研究室発． くそう， ひさしぶりに札幌駅周辺の本屋でも， と思ったけど今日は早じまいか． 2000 帰宅． 体重 73.8kg． 運動． 晩飯．
• [今日の運動]
  • エアロバイク 45 分間．
  • 腹筋運動 30 × 3 回． 腕立ふせ 5 × 3 回．
• [今日の食卓]
  • 朝 (0800): 米麦 0.7 合． キャベツ・ブナシメジ・コンニャク・油揚・煮干の味噌汁．
  • 昼 (1310): 研究室お茶部屋． 米麦 0.7 合． キャベツ・ブナシメジ・コンニャク・油揚・煮干の味噌汁．
  • 晩 (2300): 米麦 0.7 合． キャベツ・ブナシメジ・コンニャク・豆腐・油揚・煮干 の味噌雑炊．