ぎょーむ日誌 2009-05-(21-31)

2009 年 05 月 21 日 (木)

• 0830 起床． 朝飯． コーヒー． 0930 自宅発． 晴． ちょっと暑い． 0945 研究室着．
• 1000 より 1.5 時間ほど， A 棟耐震工事の疎開場所である工学部 Q/R 棟 2-3F を見にいくツアーに参加． ネット配線まわりを調べた． Q/R 棟は外観は新しいんだけど， これは最近になって耐震工事をしたためで， 内装は古くさい． そのおかげで ethernet 配線まわりは比較的にわかりやすい． ただし， R 棟 3F の hub が 2F においてあったりして， 10base-T ケイブルの限界 (50 m) をこえるケイブル長にならないか心配だ．
• ( 夕方になってから，だけど) 研究室雑用 wiki に上の調査結果などをまとめてみた． こっちにもはりつけとくか．
  • R 棟 3F R364 院生部屋
    • 現時点でネットワークなし
    • R363 または R311 とつなぐ配線を依頼 (2009-05-21)
    • 可能なマシンについては無線 LAN による接続，配線量をできるだけ少なくする
  • R 棟 3F R363 実験室
    • R312: R 棟 2F R254 室前の hub G14A9 port 08
  • R 棟 3F R312, R313A, R313B 教員部屋
    • R312: R 棟 2F R254 室前の hub G14A9 port 16
    • R313A: R 棟 2F R254 室前の hub G14A9 port 13
    • R313B: R 棟 2F R254 室前の hub G14A9 port 15
  • Q 棟 3F Q305 事務室
    • ネットワークなし，配線を依頼
    • 北大ネットの DHCP で対処する
    • LAN にはつなげないであろう
• 昼飯．
• 昨日は失敗した cairo R コンパイル作業のやりなおし． 失敗の原因は， 昨日も書いてるように pkgconfig をインストールしてなかった， なる Vine Linux に固有な問題であるように思われるので， 本日はただ単に (pkgconfig がすでにインストールされた状況で) time rpm -ba --clean ./R.spec すればよい， と． また 83 分かかってコンパイル …… cairo R ができた! X11() や png() といった device で下のような半透過 (semi-transparent) な作図ができる．

• ということで， $HOME/.Rprofile にて grDevices::X11.options(type = "cairo") その他あれこれと設定すればよい， と． しかし， Ripley 教授系の library(grDevices) な cairo device と， Simon 系の library(Cairo) が並立する現状， 将来はどうなっていくんでしょうねえ …… (後記: cairo な X11() はみょーに描画速度が遅い場合がある， とわかった． grDevices::X11.options(type = "Xlib") にしておいて， 必要なときだけ type = "cairo" にしたほうがよいかも)
• A802 院生部屋のプリンター CANON LBP5400 不調問題のつづき． 生協のヒトが実験的に M トナーカートリッジを新品に置き換えたところ …… 「M トナーカートリッジぎれ」 のニセ警告は消失した． ということで， 今まで使ってきた再生トナーカートリッジのハードウェア的な故障が 原因らしい， とわかってきた． ということで， 調査対象は本体納品した生協ではなく， 再生トナーカートリッジ屋に移っていくことになる．
• Pasoh で「いっせい」開花 (といっても 4 樹種だけ) したということで， 矢澤さんがマレイシアにぴゅーっと飛んでいってしまった．
• また北大祭関係の書類かき & 提出．
• なぜかデータ解析こんさるメイルが同時多発的に来信する． Inbox に未処理なメイルがこれ以上蓄積するのは耐えがたいので， 「その場で撃ちかえし」 な作業にはまりこんでしまう．
• それにしても random effects という概念は， なかなか理解してもらえないものですなぁ ……
• けっきょく， 今日は研究方面は何もススまなかった， とよろよろと撤退しようとしたところ， 生態学会事務局から不幸のメイルが． 和文誌メイルアカウントにきたメイルが 茨城大の編集部あてに転送されません?
• 転送設定は問題ナシ， に見える …… しょうがないので， さらに生態学会さーばー機の深奥にダイヴしてみると …… 原因がわかった． /usr/local/psa/var/log/maillog.processed の記録をみると， May 19 00:37:10 ve441 qmail: 1242661030.636380 delivery 43482: deferral: Connected_to_157.80.12.22_but_sender_was_rejected./ Remote_host_said:_450_4.7.1_Access_denied._IP_name_lookup_failed_ [203.104.99.221]/ …… えーい， つまり茨城大のメイルサーヴァーが 「203.104.99.221 (これは生態学会さーばー) とかいう DNS 逆引きができないところからきたメイルなんて受けとれるか! どうせスパムだろ!」 と決めつけてくれている， とわかった …… たしかに生態学会さーばーから送られてくるようなメイルってのは， もっぱら研究を妨害するスパムである確率が高いのはたしかなのだが ……
• まあ， 「直せ」とのことなので， 二段がまえで対処することに． くだんの IP アドレスからメイルを通すように， ひとつはその和文誌編集部が大学のネット管理者と交渉すること． もひとつは， 当方からさーばー屋さんに上記 IP アドレスを (host 名は何でもいーから) DNS 逆引きできるようにそちらさんの DNS さーばーを設定してくれよ (どうせまた似たような問題が発生するだろうし)， というもの． 後者だけでいいだろ， と思われるかもしれないけど， こっちはどれだけ時間かかるかわからず， 即効性がないのである …… 前者は編集部 vs 大学ネット管理者の抗争ぢゃなくて 交渉の進捗依存なので， 当方は無責任をきめこむことができる．
• いやはや …… 1950 研究室発． 2010 帰宅． 晩飯の準備． 晩飯．
• [今日の運動]
  • またうんどう休養日になってしまった …… このあたりちょっと発想を変えて， うんどうにとりくみやすくなるような方法を導入してみよう ……
• [今日の食卓]
  • 朝 (0840): シリアル．
  • 昼 (1230): 研究室． ハムサンドイッチ．
  • 晩 (2200): 米麦 0.8 合． ニラ卵炒飯． ダイコン・ニンジン・昆布・タラすりみの煮物． ヨーグルト．

2009 年 05 月 22 日 (金)

• 0830 起床． 朝飯． コーヒー． 0930 自宅発． 雨． 0945 研究室着．
• 院生の PC の電気鼠不調問題に 2 時間ばかり使ってしまった …… どうもこれはハードウェア的な問題ではないかな?
• 昼飯．
• またメイル書きとか．
• 昨日の生態学会さーばー DNS 逆引き問題．

  IPアドレスの所有者は、業務提携先のライブドア社となっておりますので、 現在逆引きが可能か確認をしております。

  なんと livedoor の DNS server 使ってたのか． 夕方には逆引きできるようになった．
• ZI 解説論文原稿の見なおし …… えーい， Table を新しく使って整理整頓できないかな?
• 15 時から一時間ほどこんさるぎょーむ．
• また ZI な世界にひたる． どういうふうに改訂したもんかな ……
• 1850 研究室発． 雨． 1905 帰宅． 晩飯．
• 太陽を横切るスペースシャトル (wired vision) …… うーむ， この不思議感のある「絵」をどこかで「見た」ような気がするなぁ? …… そうか! 「火星鉄道一九 (Martian Railroad 19)」 (谷甲州の航空宇宙軍史シリーズの短編集) 収録の「水星遊撃隊 (Mercury Scout)」 (1987 年) の一場面だ．

  「ちょっと待って下さい．いま観測結果が出てきました．あの艦は，やっぱり進行方向と 直角に加速しています．我々のレーダー観測時には，太陽と本艦を結ぶ線上に あの艦は入っていました．方向角の誤差は数秒くらいなもんです」
  「本当か? 偶然じゃないだろうな」
  「間違いありません．どんぴしゃりでそこにおさまっています」
  「ということは，掩蔽の観測をしやがったのか．ただ者じゃないな……．奴らは」

  (中略)

  敵は掩蔽を利用して哨戒艦の観測を行った． 光学観測は，普通なら他のセンサの補助手段としてしか使用されない． 目標が太陽を背後にしている時は， 特にそうだ． しかし， 掩蔽が起こるときは事情が異なる． 精度のいい光学望遠鏡があれば， 太陽の光球面をよぎる艦影をとらえて，速度や移動方向を知ることができる．

• [今日の運動]
  • 今日は自宅発ではなく研究室発の北大構内走しようと思って， ウェア & シューズもってきたんだけど雨だった ……
• [今日の食卓]
  • 朝 (0840): 食パン． ヨーグルト．
  • 昼 (1230): 研究室． 食パン． リンゴ．
  • 晩 (2130): スパゲッティー．

2009 年 05 月 23 日 (土)

• 0900 起床． 生活周期， ずれぎみ． 朝飯． コーヒー． 洗濯． 怠業． 1120 自宅発． 雨． 1135 研究室着．
• Atom 生成なしくみの修復， 30 分ほどじたばたしてみた．
• 来てみたものの， 怠業状態がつづく ……
• 気分転換のつもりで， 15 時ごろから一時間ほどクラーク会館の 理髪店 で散髪． 外はあいかわらず雨．
• 散髪してアタマはすっきりしたものの， アタマの中はすっきりせず．
• 現在， アタマの中にいろいろ検討すべき問題がちらかりまくっているのだが …… Pasoh 外来種の一次元空間分布問題， ベイズな空間統計モデリングをもっともお手軽にやろうとすると， car.normal() つまり intrinsic Gaussian Random Field (intrinsic GRF) で， となる． しかしながらちょっと前までの試行錯誤で， 空間 zero inflated なデータに対して， GRF のたぐいはあまり有向でないかも， という気がしてきた (上述の car.car() 例題 みたいな場合は GRF で何も問題ナシ)．
• ガウス確率場 (GRF) がダメなら …… たとえばマルコフ確率場 (Markovian Random Field; MRF) はどうだろう? MRF じたいは特に難しいものではない． ようするに イジングモデル だ． しかし観測データと対応づけるためには工夫が必要とされる． そのためには， 「誤認識確率」のようなものをうまく導入してやる必要がある， と …… 画像修復みたいな問題ですな．
• …… てなことをうだうだと考えていて， 実際には何も進捗せづ．
• 東南アジアから論文みてくださいメイルがくるし． 現在， 「論文みてください (修正せよ!)」案件が 3 つ， 放置計算したうけがひとつ (いやふたつだ)， 放置データ解析こんさるがひとつ， 放置事務処理がひとつ， 再投稿したいけど放置されてる論文原稿がひとつ， さっさと書いてしまいたい論文原稿がひとつ， さっさと修正すべき書籍原稿がひとつ， さっさととりくまねばならぬ書籍原稿がひとつ ……
• うーむ． 2045 研究室発． 雨はやんでる． 買いもの． 2110 帰宅． 晩飯の準備． 晩飯．
• [今日の運動]
  • どーも不調ですなぁ ……
• [今日の食卓]
  • 朝 (0930): 食パン． 卵焼き．
  • 昼 (1300): 研究室． 食パン． リンゴ．
  • 晩 (2140): ニラ・エノキダケ・ニンニクのスパゲッティー． ヨーグルト．

2009 年 05 月 24 日 (日)

• 0840 起床． 朝飯． コーヒー． 洗濯． 怠業． 1040 自宅発． 曇． 1100 研究室着．
• 昨日列挙してみた「放置」案件のうち， とりあえず他人の人事にかかわる放置事務書類かき …… うだうだと能率わるく進捗させたので， 2 時間ほど費やしてしまった．
• 昼飯． ひと月後ぐらいにせまった A 棟耐震工事疎開ひっこしが気になるので， お茶部屋においてた Nature バックナンバーをヒモでくくって 廃棄できるようにする． これからは毎日のようにゴミすて作業をしないといけないような気もする．
• 作業テイブルの脚がぐらついていたので， モンキーレンチでボルトをしめなおしてみる． そうか， このテイブルは解体できるのか， とまたひっこし作業について考えてしまう．
• 次． 放置データ解析である Alpinia 結実． 先週， R で自動データ点検やって入力まちがいと思われる箇所を いくつかみつけて， 発注者に「おかしいんぢゃない?」 とお尋ねしたところ， やはりまちがいだったとわかった．
• で， さっさとまちがいを修正して計算やりなおせばよかったのだが …… まあ， なンともやる気がでなくて放置していた次第．
• CSV ファイルをひねくってデータ修正 (われながらザツだなぁ)． データ変換プログラムを修正して「結実数ゼロの果実」 に対処できるようにする (じつはここでしくじって， このあと一時間ぶんの計算がムダになった)．
• ところで， 先日 cairo R を作った， といったハナシをしたけど， X11(type = "cairo") は X11(type = "Xlib") に比べてかなり描画速度が遅いような気がする． $HOME/.Rprofile は grDevices::X11.options(type = "Xlib") と設定しておいて (cairo が使えれば type は自動的にそちらを使ってしまう)， 必要なときだけ cairo したほうがよいかも．
• うーむ， library(Cairo) の CairoX11() は多少速い (?) ような気がする．
• で， 個体差と個花差が 「ホントに個体群なのか?」 と思えるほどに絶望的にばらばらなデータなので， 絶望的に計算時間が長くなってしまう階層ベイズモデルの推定やりなおし．

[個体差・個花差]

あとからみつけやすい場所に図メモ． 全データの PDF ファイル．

• で， 絶望的に長い計算時間はとっとと別のぎょーむに取りくめばよいのだが …… またぎょーむ日誌の atom.xml を改良するとゆーどうでもよいプログラミングにはまってしまった． 本日は各えんとりーの <updated>2009-05-24T16:26:40+09:00</updated> タグを改善する， といった案件で …… 自分で構築したシステム (ヒトことで言えば，テキストベイスのデータベイス) なんだけど， このあたり不必要に複雑ぎみなので 改良作業にも時間がかかる．
• Alpinia 再計算の結果でた． データ修正個所が多くなかったせいか， 推定結果 は 以前 とほぼ同じ， と．
• 1830 研究室発． ちょっとうろうろ散歩． 晩飯． 2030 帰宅．
• Alpinia したうけ， 一段落したと思ったら論文かけ指示がでてがっくり． どんどん「論文修正しろ」「論文かけ」 が堆積していく ……
• [今日の運動]
  • なかなかうんどう再開できない日々がつづく ……
• [今日の食卓]
  • 朝 (0920): 食パン． ニンジン・ニラ・卵の炒めもの．
  • 昼 (1330): 研究室お茶部屋． 食パン． ワカメスープ． リンゴ．
  • 晩 (1930): 心 のスープカレー．

2009 年 05 月 25 日 (月)

• 0810 起床． 朝飯． コーヒー． 0940 自宅発． 晴． 1000 研究室着．
• 評価書なる事務書類の提出． 6 部コピーをとって角 2 封筒にいれて密封して， 学術助成係に提出， と．
• 夕方までひたすら ZI 解説論文の修正作業． MS わーど (on VMWare) を使わねばならないので， 不幸なキモチになる． わーどを使わなくてよい論文かき (すなわち自分が first author) に復帰できるのはいつの日になるのだろうか ……
• とりあえずひととーり修正できた． 沓掛さんに送信 …… 時刻はすでに 18 時すぎか．
• 院生たちの圧力を感じたので， お茶部屋で「捨てるべき本」の選別を少しばかり． これからひと月， ゴミ分類の日々がつづきそう ……
• 「うんどう休養日」つづきに対処すべく， 自宅発ならぬ研究室発北大構内走を実施してみる． とにかく帰宅するとぐったりしてしまって， 何もする気になれないんだよね …… 1845 研究室発北大構内走． ちょうど日没時刻． 1925 研究室もどる． けっこう汗がでてきた．
• さっさと帰ることに． 1940 研究室発． 買いもの． 1955 帰宅． 体重 72.8 kg． 晩飯の準備． 晩飯．
• 電気代 1439 円． 毎月のばらつきは何に依存するのだろうか?
• [今日の運動]
  • 北大構内走 40 分間． ストレッチ．
• [今日の食卓]
  • 朝 (0850): 食パン． ヨーグルト．
  • 昼 (1250): 研究室． 食パン． リンゴ．
  • 晩 (2130): 米麦 0.8 合． キャベツ・ニンジン・エノキダケ・鶏肉の煮物． ヨーグルト．

2009 年 05 月 26 日 (火)

• 0600 起床． 朝飯． コーヒー． 0725 自宅発． 晴． 0740 研究室着．
• セミナーの予習．
• 1030 より 研究室セミナー， 本日は川久保さんで火山噴出物・種子散布・窒素固定植物・地温の 4 論文の紹介． このうち 3 つは Journal of Vegetation Science なのだが …… 空間構造・時系列・多種系というもっとも複雑なる生態学データを とっておきながら， 可能なかぎりヘンてこにデータをあつかう論文を掲載するぢゃーなる …… とゆー気がしてきた． 火山とか極端環境でうまい研究をするのは， 案外むずかしいのではないかな．
• 昼飯．
• 一次元空間モデリング， 試行錯誤のつづき． 脱・Gaussian random field (GRF) のため有限状態の Markovian Random Field (MRF) を検討してみる (注: 時系列データ解析方面のモデルと混同しそうでややこしいのだが …… MRF の定義は「各地点の状態数が有限」ではなく， ある地点に影響を与える「近傍集団」が有限の範囲におさまる， ということ． したがって WinBUGS の car.normal() な GRF もまた MRF の一種である)．
• 出発点としては， 各地点において {0, 1} の二値をとる MRF を考える． これを観測データにあわせるためには， 「観測をまちがえる確率」 みたいなものを導入しなければならない …… そしてこれは観測データから推定することは不可能， とわかってきたので P.err = 0.001 などと決めてやらねばならない ……
• で， とうぜんながらというか， 観測データと対応させようとすると， 各地点の値はほとんど動かなくなり (遷移は上述の P.err にほぼ支配されている)， あんまし確率場 (random field) ってかんぢではなくなる ……
• そこで， べるぬーい (dbern(p[i])) からばいのみある (dbbin(p[i], N.size)) に拡張してみる ……
• こういう新しい試みはなかなかうまくいかない …… とりあえず， 途中で一時間ほど藤沼君と現状について検討してみる． 他人と検討してみると， 少しアタマが整理される．
• また別のデータ解析こんさる 40 分間ほど．
• 来週の土日の北大祭雑用， 静かなる混乱というべき状況にあるようだ． 担当者と事務方面 (全体ならびに各専攻) がばらばらに動いている． 全体の事務から必要なら出勤とどけをだせ， というメイルがきたんだけど， とくに仕事は assign されていないので， そういう紙きれは出さないことにする． まあ， どうせ土日は大学にいるだろうけど …… 土曜日は勉強会もあるし．
• 二項分布な MRF， 試行錯誤しながらその性質を理解しようとしている …… だんだんわかってきた．
• 1845 研究室発北大構内走． てけてけ走っているとアタマの中が整理されてくる …… 1930 研究室もどる．
• ということで， これでうまくいくだろう， と思いきってパラメーター数を減らしてみたら， うまくいった! 翌日以降， またいろいろと変更するかもしれないので， とりあえずここに記録しとこう (BUGS code， 事後分布表)．

• さて， この程度のことが GRF ではうまくできないのか， と再確認してみたくなった． まず Distance 依存性がない最単純モデルによる推定 (BUGS code)． この程度の計算は， さすがに収束する．

• ちょっと横道にそれるが， 私の場合， こういう単純なモデルによって GRF のような (じつはけっこう難しい) モノの性質の理解が深まる． たとえば上の binomial なやつとこの GRF なモデルの相違は， intcpt なる「切片」パラメーターの有無にある． GRF なほうでは intcpt ナシではうまく収束せず， binomial MRF なやつだとこういう「切片」パラメーターは ただのジャマもの (ということに北大構内走の最中に気づいた)． intrinsic Gaussian CAR normal なモデルは計算はカンタンなんだけど (条件つき局所確率は単純な正規分布にしたがうから)， その「実態」の理解は私などにはなかなか難しい (確率場全体を簡単な多変量確率分布で定義するのが難しそう) といったシロモノなのだが …… もしこの car.normal() で生成される GRF をひとつの多変量正規分布 (?) で書けるとすると， その「各地点の平均はことごとくゼロ」 と設定されているはずだ， とわかった (ただし分散共分散行列のほうはナゾで， 定義できないのではないかと思う …… まあこのあたり variogram と対応可能な GRF あたりから勉強しなおすべきかも)．
• で， 上の単純な GRF モデルで気になるのは， 気になるのは「zero inflated な地帯」ともうしますか， ひたすら -0-0-0-0-0-……と続いてる場所において (ひたすら -1-1-1-1-1-…… でも同様だろうが)， 状態変数が -10 から 0 ぐらいの範囲をぶんぶん 振動しまくっていること． これはまったく妥当な事後分布なのだが ……
• しかし， このあたりが原因 (?) で， モデルをちょっと複雑にしたときに収束しなくなるのかも， という気もする． 今度は Distance 依存性も推定させようとしたもの (BUGS code， 事後分布表)． 上の binomial なやつではまがりなりにも距離依存性を 推定できていたわけだが， こちらの GRF なやつだとかくのごとくコケてしまいます， と．

• 上のように計算しながら， 院生のデータ解析コンサル． R による作図と glm() によるやや複雑なモデリング …… 「あろめとりっく」なハナシなので， 説明変数を log.x としながら， family = gaussian(link = "log") とするわざ (今年 3 月の盛岡での GLM 自由集会 の久保投影資料を参照)．
• もちろん， こんなごまかしわざ (盛岡での粕谷さんのポスター発表で この「x と y に『誤差』がある場合の回帰の問題点」 が指摘されてましたね) ではなく階層ベイズモデル化してケリをつけたいところだが …… まあ， 院生たちはまだゐんばぐすは修行中， ということで． しかし， 推定結果と観測データを比較してみたら， やっぱり「個体差」がすごかった．
• 2145 研究室発． 2200 帰宅． 晩飯の準備． 晩飯．
• [今日の運動]
  • 北大構内走 45 分間． ラジオ体操． ストレッチ．
• [今日の食卓]
  • 朝 (0700): 米麦 0.6 合． キャベツ・ネギのサラダ．
  • 昼 (1230): 研究室． 食パン． チーズ． キャベツ・ネギのサラダ．
  • 晩 (2230): 米麦 0.8 合． キャベツ・ニンジン・エノキダケ・鶏肉の煮物． サバの電磁波酒蒸し． ヨーグルト．

2009 年 05 月 27 日 (水)

• 0840 起床 …… うーむ， これは昨日の労働時間が 14 時間超であったため， ということいしておこう． 朝飯． コーヒー． 0940 自宅発． 曇． 0955 研究室着．
• で， じつは昨日理解できたあれこれを本日になって ぎょーむ日誌にまとめていたりしている．
• 前後とまったく関係ないハナシだけど， REML method (restricted maximum likelihood method; 制限つき最尤推定法) のもとでのモデル選択ってまだまだケリがついてない ハナシ のようだ． そして， (何でもかんでも正規分布ならともかく) random effects のあるモデル (つまり mixed model) 状況下での分散推定は 難しい ……
• で， 一次元空間上での MCMC ダンス 試行錯誤のつづき． やはりゐんばぐす car.normal() では説明変数的なパラメーター数の個数を減らしてもダメ， とわかった．

• ねんのために書いておくと， ダメなのは「空間構造あり，かつ zero inflated なデータのモデリング」 だけなのであって， この例題 のようにゼロが多くないデータならばとくに問題はなかったりする．
• 雑用がときどきわりこむ． 研究室ネット雑用． 研究室ひっこし雑用の相談． などなど．
• そうそう， 北大にしばらくいる Paul Moorcraft のプリンター設定もやった． 会うのは 10 年以上ぶりですなぁ ……
• auto-logistic (binomial 版) の試行錯誤つづく． この新しいモデルでは， こうしてみるのが妥当なんだろうなあ， というのがわかってきた (BUGS code， 事後分布表)． dbin() の size を増やすと (みかけのうえで) ばらつきは減る． ただし計算には時間がかかる (二項乱数の生成を考えるとそうかも)． 1200 秒ぐらい．

• などといろいろやってるうちに， 一次元空間上の確率場な統計モデルで 「あーすればこーなる」 というのがわかってきた．
• ということで， ここ数日のあいだに， 何度も見捨てられては見なおされてきた Gaussian Random Field (GRF) を使ったモデリングの再検討 …… しているうちに， 私の耳もとで誰かが 「あれでしょ? よーするに， こーすれば機嫌よく収束してくれるんでしょ?」 とささやいて， (ぢつに直観的な) 小細工をくみこんだ新しい定式化を ぽいっと渡してくれた …… 使ってみたらうまくいった． GRF なのにすぱっと収束， しかも計算時間はわずか 50 秒に短縮された．

• えーと， やりかたとしては， intrinsic Gaussian CAR だからといって， できあいの car.normal() 関数にたよるのではなく， BUGS code 内で自分で書いてしまう， ということ． それによって小細工を混入させうる余地が生じる． 小細工のカギは局所的な条件つき確率の定義にあり，

  	x[i] ~ dnorm(q[i], tau)
  	q[i] <- (x[LeftRight[i, 1]] + x[LeftRight[i, 2]]) * 0.45
  

  の 2 行目で 0.5 ではなく 0.45 にしている， というところ． 0.5 にすると car.normal() と同じ動作になってしまう (そして収束しない)．
• さてさて …… このように「平均」をとるのが合理的だっていうのは， 上の推定計算結果を見てしまったあとでは， 「まあ，そうでしょ」となるわけだが …… 数学的にはどう正当化されうるのかよくわからない．
• いやはや， 私の場合， イチローのように 「持っている」 わけではないので …… 神様ならぬ統計モデリングの気まぐれな悪魔が渡してよこした， この格段に良い方法候補についてよくよく検討してみる必要がある．
• これがホントに正しい方法であるなら， 誰かがすでに理論的な検討をしていそうな気もするのだがなぁ．
• よくわからないので， 走ってみることに． 1945 研究室発北大構内走． 「まあ， あのように計算すればうまくいくよなぁ」 といった理由はいろいろ思いつくんだけど …… 「うまくいくから，よくわからんけど OKOK」 というのはかなり危ない考えかただ． 少なくとも理学部な世界では． もうちょっと原理・原則にもとづいたうまい導出が …… 今のところ思いつかない． 2030 研究室もどる．
• 2050 研究室発． 2115 帰宅． 体重 73.0 kg． 体重はなかなか減少しないけど， まあ気長に． 晩飯の準備． 晩飯．
• [今日の運動]
  • 北大構内走 45 分間． ストレッチ． さいきん， カラダがカタいような気がする．
• [今日の食卓]
  • 朝 (0910): 米麦 0.6 合． サバの電磁波酒蒸し．
  • 昼 (1220): 研究室． 食パン． チーズ．
  • 晩 (2230): 米麦 0.8 合． ダイコン・ニンジン・ブナシメジ・鶏肉のコンソメ煮物． ヨーグルト．

2009 年 05 月 28 日 (木)

• 0700 起床． 朝飯． コーヒー． 洗濯． 一次元空間上の確率場についてうだうだと考えてみる． 0930 自宅発． 晴． 川むかいの北 12 条郵便局に転送してもらった 定額給付金の紙束をうけとる (札幌市は全国の自治体の中でこの紙束を郵送するのが もっとも遅い)． まあ， 15 兆円の補正予算の一部はのちのちわれわれが支払わされるので (納税者ひとりあたり 20-30 万円ぐらい?)， 1.2 万円ぐらいは返してもらってもいいでしょう， ということで． 0955 研究室着．
• 昨日のつづきで， また一次元 GRF 問題で試行錯誤してみる． しかしあまり進捗なく，

  	x[i] ~ dnorm(q[i], tau)
  	q[i] <- (x[LeftRight[i, 1]] + x[LeftRight[i, 2]]) * w
  

  としたときに
  • w が 0.5 より小さい正の値であれば (そんなにいろいろ調べたわけではないが)， まあだいたい同じような結果に収束する
  • w = 0.49 など 0.5 に近づけていくと， 徐々に収束が遅くなる
  • w = 0.5 は「臨界点」みたいな値で， この値だと一気に収束が遅くなる (あるいは永久に収束しない?)
• 雪野さんとプリンター相談， とか．
• で， 夜まで藤沼君と相談しながら Pasoh ブタ-Clidemia 空間統計モデリング． 段階的に試行錯誤しながら， いろいろな要因をほうりこんでみる．
• とりあえず， ひととーりあれこれを考慮したモデリングができるようになった …… とはいえ， これにはまだまだ改善すべき点がある， とあとになって気づいた． (BUGS code， 事後分布表)． 計算時間は 1300 秒ぐらい．

• 1915 研究室発北大構内走． こんな統計モデリングでよいのだろうか， と検討しつつ． 1955 研究室もどる． ばてた． やはり他人にあれこれ説明しながら， 試行錯誤していたのでアタマがクタびれてた．
• 2030 研究室発． 買いもの． 2050 帰宅． 晩飯の準備． 晩飯．
• [今日の運動]
  • 北大構内走 40 分間． ストレッチ．
• [今日の食卓]
  • 朝 (0730): 米麦 0.6 合． ニラ卵炒飯．
  • 昼 (1240): 研究室お茶部屋． 食パン． ダイコン・ニンジン・ブナシメジ・鶏肉のコンソメ煮物．
  • 晩 (2200): 米麦 1.0 合． ダイコン・ニンジン・豚モツのカレー． ヨーグルト．

2009 年 05 月 29 日 (金)

• 0830 起床． 朝飯． コーヒー． 0940 自宅発． 晴． 0955 研究室着．
• また Pasoh ブタ-Clidemia 一次元空間統計モデリングの改良のつづき． とりあえず「説明変数」(といっていいのか?) を中央化 (centralization) したら計算時間が半減して 700 秒になった．
• 昨日・本日， いろいろな定式化を考案しては試行錯誤しているのだが …… 空間統計モデリングの細かいワザ・要点あれこれがわかってきた． たとえば， fixed/random effects と分類する世界とはちょっと違うのだけど， たとえば生物の密度予測の問題において， 線形予測子の構成を絶対値・相対値に分割， つまり
  1. 「切片」 (intercept) 的なパラメーター (例: 調査地全体の対数平均密度の値)
  2. GRF で生成される値 (この GRF は説明要因によって局所平均値が上下する)
  にわけて考えるのが良いような気がする． で， ポイントはいわゆる説明変数的な因子も GRF 生成のほうに ほうりこむのが正しいように思える． なぜかというと， たとえば「ある地点だけ明るい」→「その地点だけ生物密度が明るくなる」 →「局所的な移動によってその周辺の密度も少し高くなる (空間的自己相関の発生)」 となるからだ． fixed effects は線形予測子で， random effects は GRF でといったわけかたをしていると， このあたりがうまく定式化できない．
• もちろん局所移動とかが考えられない生物で， 生物の「移動」による空間的自己相関はまったくない状況では， GRF 計算の中に説明変数をいれるのはまずいよ．
• もひとつ補足しとくと， これは空間上の要因を統計モデル化する必要がない， あるいはデータが少なすぎてモデル化不可能， という場合の措置だ． このブタ-Clidemia 系だと Gap だの Swamp だのがそれに該当する． あまりにも個数が少ないのでモデル化はほとんど無理． いっぽうブタによる撹乱は明示的に空間統計モデル化しており， 撹乱の集中度や環境要因の影響が推定される， と．
• またべつの工夫． これは空間統計モデリング限定ではなく， 複数の状態変数を同時に推定する統計モデリング全般にいえることなんだけど …… たとえば，

  	Nobs[i] ~ dbin(p[i], N.size)
  	logit(p[i]) <- a + xN[i]
  	Mobs[i] ~ dpois(q[i], N.size)
  	log(q[i]) <- b0 + b1 * xN[i] + xM[i]
  

  といった BUGS コードで，
  • i: 場所
  • Nobs[i], Mobs[i]: 観測値
  • xN[i], xM[i]: 潜在変数 (状態変数)
  となっているときに，

  	log(q[i]) <- b0 + b1 * xN[i] + xM[i]
  

  といったかんじで， 「状態変数 xN[i] をそのまま 状態変数 xM[i] の説明変数に使う」 というのはマズい場合があるかもしれない． xN[i] が xM[i] のあてはめに 「ふりまわされる」 ような状況になることがある．
• このような状態変数間の相互作用を弱めたい場合には，

  	Nobs[i] ~ dbin(p[i], N.size)
  	n.simulated[i] ~ dbin(p[i], N.size)
  	logit(p[i]) <- a + xN[i]
  	Mobs[i] ~ dpois(q[i], N.size)
  	log(q[i]) <- b0 + b1 * n.simulated[i] + xM[i]
  

  このように n.simulated[i] なる補助的な状態変数 (N なる量を ``simulate'' した量) を説明変数にすると， xN[i] が xM[i] にあまり影響されなくなる． 影響されないほうが良い場合には， このように定式化したほうがよい． (そして n.simulated[i] も「それっぽい」値を引いて「なんちゃって中央化」 をすると収束が改善される)．
• 本日また確認できたこと． 一昨日に比べてかなり複雑になった一次元空間モデルなのだが …… 依然として，
  • GRF の局所確率分布の平均を (x_left + x_right) * 0.5 で評価すると収束がむちゃくちゃに遅くなる (もしくは永久に収束しない?)
  • たとえば (x_left + x_right) * 0.45 といったナゾのやり口だと問題なく収束しているように見える
  といった状況にある， と確認できた．
• 計算まち時間には雑用とか．
• 来週の研究院教授会議の案内がきたので， そくざに「欠席」の返信．
• で， 上の要点あれこれを反映させると， まあそれっぽい結果が得られました， と (BUGS code， 事後分布表)． 計算時間は 700 秒ぐらい．

• で， 藤沼君あいてにモデリングの方法と結果についての説明． とりあえず， 上の結果にもとづいて論文原稿を書いていくことにする．
• ちょっとデータ解析こんさるメイルかき， とか．
• 5/14 以来， かなり苦闘してきた空間統計モデリングがひとまずは終了したせいなのか …… やや呆然としている．
• 呆然としている場合ではないので， 現状を少し整理してみよう．
  • 宮田さん苫小牧樹木シュート論文かき: 現時点は宮田さんターンなのだが， ご本人が耐震工事疎開ひっこし・GCOE その他雑用ストーム にふりまわされているので， なかなか進捗してないようだ
  • 矢澤さん熱帯林樹木カタチ論文かき: 当方に「原稿みてくれ」と送りつけられた状態なんだけど， ご本人はまだ赤道近辺にいる． 原稿ざっと見たんだけど， 計算ぜんぶやりなおしたほうがいいような気がするんだけどなぁ …… てなコトを言ってみたらどうなってしまうのかしらん?
  • 藤沼君一次元空間モデリング論文かき: まあ， これは上述の結果をもとに藤沼君ターン， ということで．
  • ZI 解説論文かき: 沓掛さんターン．
  • R Graphics 本の翻訳原稿の修正: 久保作業は放置されたまま …… もう 2 週間 が経過してしまったか． さっさとできる作業はやってしまおう． そして共立出版の LaTeX マクロプラミングわざ で「さくいん」が作れるかどうかが調査されているはず， なのだが ……
  • 岩波本: さっさととりかかりたい!
  • 久保 first author 論文: さっさとアリ論文を再投稿して「脱」アリしたい， さっさとウミガメ論文を書いてしまって「脱」カメしたい， さっさと ……
  • その他，放置データ解析こんさる (メイル) 一件， 放置雑用いくつか ……
  ふーむ ……
• とりあえず， 放置してたデータ解析こんさるメイルかき二件 …… これで放置案件の個数が 2 個へった．
• 1500 から面談データ解析こんさる， 50 分間ほど． 極値分布 理論に関するもの． これはなかなかやっかいだ．
• 1600 から理学部 5 号館 913 号室で， 松島さん ・川森さんヒヨコ師弟とヒヨコうろうろデータ解析の雑談． 動物認知脳科学の研究室でも GLM やベイズといった統計モデリングの方法が 普及してきたとのことで， めでたいことです． 話してるうちに， だんだん川森さんモデルの内容を思い出してきた ……
• 研究室にもどって， ちょっと空間統計モデリングの勉強．
• 1905 研究室発． 1920 帰宅． ちょっとうんどう． 晩飯の準備． 晩飯．
• The R Journal がでた． これは以前から続いていた R News が発展したもの． たしかにペイジ数とかふえてる．
• [今日の運動]
  • 腹筋運動 30 × 3 回． 背筋運動 30 × 3 回． 腕立ふせ 10 × 3 回． ストレッチ．
• [今日の食卓]
  • 朝 (0850): 米麦 0.6 合． ダイコン・ニンジン・豚モツのカレー．
  • 昼 (1220): 研究室． 食パン． 野菜ジュース．
  • 晩 (2200): スパゲッティー． ブナシメジ・ニンニク・ホタテのトマト煮． ヨーグルト．

2009 年 05 月 30 日 (土)

• 0800 起床． 朝飯． コーヒー．
• 午前中はぐったり怠業でオワってしまった ……
• 午後はだらだらと掃除・料理・洗濯 ……
• 雨． うーむ ……
• [今日の運動]
  • うんどう休養日．
• [今日の食卓]
  • 朝 (0900): 食パン． ブナシメジ・ニンニク・ホタテのトマト煮．
  • 昼 (1230): イソップベイカリーのバゲットなど．
  • 晩 (1920): 米麦 0.8 合． レタス・ネギ・ショウガ・アサリの炒めもの． ゴボウのバルサミコ炒め． ヨーグルト．

2009 年 05 月 31 日 (日)

• 0830 起床． 朝飯． コーヒー． 洗濯． 怠業．
• 1100 自宅発． 雨． 1130 研究室着．
• ちょっと銀行の ATM まで． 往復一時間ぐらいの雨のなか札幌散歩． 空間統計モデリングについて検討しつつ …… ゐんばぐすの car.normal() とかで生成される intrinsic な CAR normal な確率場， 大域的な Gibbs 分布はごちゃごちゃしてるけど， 局所的な条件付き確率分布は簡単であり …… とナニやら神秘的な確率モデルってかんぢですが， よくよく考えてみると …… 疑似尤度 (pseudo likelihood のほう，quasi ではなく) 的なとらえかたでいくと， 何のことはない
  • 平均値が線形予測子であたえられ (たいていの場合，link = identity)， さらに「近傍」にも必ず影響される (このあたりが pseudo likelihood 的)
  • 各地点「独立」な random effects をあらわすホワイトノイズ
  といった疑似尤度な GLMM と考えることができる． 疑似尤度なので最尤推定は事実上不可能であり， Gibbs 分布を MCMC することでケリをつけます， と．
• まあ， こう考えても何かが変わるわけでもないのだが， 「ならば，このあたりは好きなようにヒネくっても問題ないな」 といった確信は得られる． とにかくベイズなモデルについて検討してみるときは， 「この統計モデルの Gibbs 分布は何なのか?」 を考えてみると， すっきりする場合もある， と …… 研究室にもどって昼飯．
• ということで， 一昨日でケリをつけた (つもりになっていた) Pasoh ブタ-Clidemia 一次元空間統計モデリングをちょっとひねくってみる …… 例の「左右の平均をとる」でやっていたごまかしを， 事前分布の導入というごまかしに置喚する試行錯誤．
• 日曜日でもデータ解析こんさるはわりこむ． 個体差・(実験用) シャーレ差がないときは glmmML() の最尤推定は収束しない (sigma が無限に小さくなるから)． glmmML 紹介 にも追記．
• 計算まち時間は苫小牧シュート論文の修正作業．
• 一次元空間統計モデリングの改良は， 「あと一歩」の状態で終わった． 明日ぐらいにはなんとかなるだろう． 2020 研究室発． 雨． 買いもの． 2045 帰宅． 晩飯の準備． 晩飯．
• [今日の運動]
  • うんどう休養日になってしまった ……
• [今日の食卓]
  • 朝 (0900): 食パン． ヨーグルト．
  • 昼 (1330): 研究室． 食パン．
  • 晩 (2130): 米麦 0.8 合． ニラ麻婆豆腐． コマツナごまあえ． ヨーグルト．