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

2006 年 12 月 21 日 (木)

• 0630 起床． 朝飯． コーヒー． 0810 自宅発． 晴． 0825 研究室着．
• 昨晩帰宅する前に命じておいた苫小牧 13 樹種シュート構造の 階層ベイズモデル 60000 MCMC step 計算， 6500 秒ほど費して終了していた． モデルもずいぶんとマトモになったおかげなのか， これだけ計算させるとさすがによく収束している．

• bugs() 関数の返してきた bugs オブジェクトの summary() はこんなかんぢで． 樹種 id 番号は下記のごとし， と．

  > rbind(n = summary(d$Sp), spc.id = 1:length(levels(d$Sp)))
         Am Ap Bp Cc Fl Mk Mo Oj Ps Qm Sa Tj Ud
  n      41 52 48 34 37 35 33 45 29 37 41 37 24
  spc.id  1  2  3  4  5  6  7  8  9 10 11 12 13
  

  おもしろそうなのはこのあたり (表) か? いや， 「表」では味気ないから nested 事後分布図にしてみるか．

• これは 去年の今ごろ 苦闘していた牛原さんの屋久島ベイズモデリングの library(grid) 作図プログラムの流用なんだけど， こいつの改造に午前中ぜんぶを使ってしまった …… いや， 小山さんのサロベツ種子トラップ解析計算もちょっとお手伝いしました．
• 昨日とのちがい …… がある． つまり重量分配問題に葉重量→葉面積問題を追加すると， 重量分配問題の結果が影響をうける， といういっけん不思議な関係で． これは葉面積というデータを「追加」したことで， 葉重量モデルが補正された， と考えればよいのだろう．
  • (結果というより現象そのものが) ややアヤしげな気がしてきたんだけど， 重量分配の樹高依存性は 依存・非依存樹種に分割された (図をみるとそうかもという気はする)
  • シュートサイズ依存性はかわらず， か
  • 新しく追加された葉重量→葉面積の樹高依存性， これってただ単に 「林道わきとかでも林冠上ほど明るいわけではない」 ってことなのでは ……? (疑惑)
  • これが種差ナシなのはおもしろいけど
• 授業おわって帰ってきた宮田さんに結果の説明． 講評 「この結果で生態学会松山大会の要旨 (1/9 締め切り) が書けそうなので， とりあえず計算はこれでよろしい」 ということで， 年内の苫小牧シュート構造階層ベイズモデルの MCMC 計算からはめでたく釈放された! そう， この突貫工事的 / 敵中央正面突破的 MCMC 計算はかかる地環研 A 棟 8F 独特の政治的圧力のもとで 強行邁進させられてきた計算ぷろぢぇくとだったのである ……
  「ところで要旨にはベイズ推定やったとか書くんですか」
  「書かないとだめでしょう」
  「まだベイズ推定とかよくわかってないんですが」
  「だいじょーぶですよ． 今回の大会に提出される要旨の中で 『分散分析やりました』とか書いてるヒトたちの過半数は 分散分析の何たるかをほとんど理解してませんから」
  「と言われるとじつは私の卒研なんかでも思いあたるフシが」
  「ということで， 今回はわれわれ二人して 『ベイズがわかったふり』をよそおいつつ要旨を書いてしまいましょう」
• まだまだモデリング・計算方法などに改善すべき点が 多々あるように思えるけれど …… ともあれ終わった終わった， ということで …… といったん帰宅して洗濯しつつ昼飯． 1340 研究室もどる． 宮田さんに計算結果や R & WinBUGS コードなど引きわたして一件落着．
• 矢澤さんと Pasoh 熱帯林雑談． うーむ， Pasoh データ管理体制って ……
• 12/25 (月) しめきりの査読原稿よみにとりかかる． これが今年最後の査読， のはず．
• あまり進捗せず撤退． 明日中には終わらせたい． 1900 研究室発． 1920 帰宅．
• [今日の運動]
  • 腹筋運動 30 × 3 回． 腕立ふせ 5 × 3 回．
• [今日の食卓]
  • 朝 (0700): 米麦 0.7 合． ダイコン・ネギ・豆腐の味噌汁．
  • 昼 (1310): うどん． ダイコン・ネギ・豆腐の味噌汁．
  • 晩 (2030): 米麦 0.7 合． レトルトパウチド八宝菜． 海藻スープ．

2006 年 12 月 22 日 (金)

• 0750 起床． 朝飯． コーヒー． 0930 自宅発． 曇． 0935 研究室着．
• 査読報告かきのつづき．
• 志水さんと地衣類相談． ちょっと前の 計算 なんだけど， われながら記憶がぼろぼろに脱落している， というか． 思えば 10 月後半のこの計算いらい， ひたすら WinBUGS をあたかも「実戦化のための訓練」 であるかのよーにひたすら走らせ続けてきたなぁ ……
• 今日はいろいろな用件とかで院生密度が低く， 静寂だ．
• わりといいかげんな査読報告を書いてしまう． 1245 送信． 一件落着． 昼飯． 人口密度が低いのでのんびりと．
• R-help のやりとり見ていて思い出したんだが …… R の線形混合モデル・一般化線形混合モデルを計算してくれる library(lme4) の lmer() の結果はそれぞれ lmer・glmer クラスのオブジェクトになるのだけど， これらをベイズモデルとみなして MCMC 計算やってくれる mcmcsamp() という関数も library(lme4) に同梱されている．
• たとえば

  fit1 <- lmer(y ~ (1 | group), family = binomial, method = "Laplace", ...) # glmer オブジェクト
  mcmc1 <- mcmcsamp(fit1, n = 500) # library(coda) の mcmc オブジェクト
  

  で何やら MCMC 計算やってくれるんだけど …… 実験してみると random effects の事後分布の計算， 何やってるのかよくわからんなぁ．
• help(mcmcsamp) によると，

  The prior on the fixed effects parameters is taken to be locally uniform.

  fixed effects は何やら一様な (局所的) 無情報事前分布で， random effects は

  The prior on the variance-covariance matrices of the random effects is taken to be the locally non-informative prior described in Box and Tiao (1973). Conditional on the current values of the random effects these are sampled from a Wishart distribution.

  平均ゼロの多変量正規分布から得られている， と． その分散共分散行列の無情報事前分布 (つまり超事前分布 hyper prior) は逆 Wishart 分布． まあ， ふつーの階層ベイズモデルっぽいことやってるわけで．
• えーと， てきとうに例題を作って

  fit1 <- lmer(y ~ (1 | group), family = binomial, method = "Laplace", ...) # glmer オブジェクト
  

  をやると

  Generalized linear mixed model fit using Laplace 
  Formula: y ~ (1 | group) 
   Family: binomial(logit link)
    AIC  BIC logLik deviance
   1385 1395   -690     1381
  Random effects:
   Groups Name        Variance Std.Dev.
   group  (Intercept) 0.0353   0.188   
  number of obs: 1000, groups: group, 10
  
  Estimated scale (compare to  1 )  0.9976 
  
  Fixed effects:
              Estimate Std. Error z value Pr(>|z|)
  (Intercept)    0.108      0.087    1.24     0.22
  

  となり， さらに

  mcmc1 <- mcmcsamp(fit1, n = 500) # library(coda) の mcmc オブジェクト
  

  やってから (library(coda) も忘れず実行) summary(mcmc1) やると

  Iterations = 1:500
  Thinning interval = 1 
  Number of chains = 1 
  Sample size per chain = 500 
  
  1. Empirical mean and standard deviation for each variable,
     plus standard error of the mean:
  
                     Mean     SD Naive SE Time-series SE
  (Intercept)       0.109 0.0999  0.00447        0.00528
  log(grop.(In))   -3.391 1.1070  0.04951        0.14390
  deviance       1385.842 6.2637  0.28012        0.47325
  
  2. Quantiles for each variable:
  
                     2.5%       25%      50%      75%    97.5%
  (Intercept)      -0.084    0.0474    0.104    0.169    0.318
  log(grop.(In))   -6.022   -3.9808   -3.269   -2.629   -1.521
  deviance       1375.724 1381.4101 1385.200 1389.406 1399.566
  

  で fit1 と summary(mcmc1) を比べてよーやくわかったんだけど， MCMC 計算の結果に含まれてる log(grop.(In)) ってのは group の random effects (つまり group 差) あらわす事後分布の分散の log になってるわけね． plot(mcmc1) すると， いつものごときこういう図になるわけで．

• つまり線形混合モデル・一般化線形混合モデルの階層ベイズモデルの MCMC 計算は R2WinBUGS なんぞ使わずとも lmer() & mcmcsamp() 使えばできるらしい， ということですな． ふーむ …… とやや驚いたので 「生態学のデータ解析 - lmer MCMC 計算」 なるペイジを追加．
• で， ついでについついこのデータ解析サイトのあちこちも変更してしまう．
• 私はなぜかくもうだうだしてるのだろうか?
  1. 今年最後の査読ぎょーむが終了したから
  2. 院生密度が低くてのんびりできるから (年内に完了すべき院生したうけぎょーむはまだ残存してるが)
  3. 12/27 (水) しめきりの 助教書類 書くのがイヤだから
• とも言ってられないので書類の準備． 途中でイヤになってお茶部屋に逃げようとしたら， 院生密度がいつのまにか平常の値にもどっていたので， こそこそと逃げだす．
• 完成度 40% ぐらい? で本日は撤退． 1920 研究室発． なんかホントに雪がぜんぜんつもってないな …… 1935 帰宅． 体重 70.6kg． 晩飯の準備． 晩飯．
• この ひと月 でさらに体重 1kg ほど減少したのはいいんだけど (下うけやせ?) …… どうも空腹状態で晩飯くうと カラダの血液のたいはんが胃のまわりにあつまるかんぢですなぁ． アタマが働かない．
• [今日の運動]
  • 腹筋運動 30 × 3 回． 腕立ふせ 5 × 3 回．
• [今日の食卓]
  • 朝 (0820): 米麦 0.5 合． 海藻スープ．
  • 昼 (1300): 研究室お茶部屋． 食パン．
  • 晩 (2110): 米麦 0.7 合． マイタケ・ワカメ・ホタテ・豆腐のシチュー．

2006 年 12 月 23 日 (土)

• 0730 起床． 朝飯． コーヒー．
• 昨晩のことなんだけど， 生態学会の 松山大会 の大会企画委員会 (委員の皆様，ホントにおつかれさまです) からシンポジウム・フォーラム・自由集会一覧が送られてきた (シンポジウムなど企画者あてに)． これってまだ公表用のものではないんだろうけど， あればお役にたつヒトもいるだろうから， このあたりにはってみる． あくまでも， 参考程度に， ということで．

● シンポジウム・フォーラム・自由集会 一覧

◎ 3/19 (MON)
自由集会 (18:00-20:00) [Room]
W01 生物群集と生態系をつなぐ：群集生態学の新たなる挑戦 [A]
W02 外来植物の侵入実態と生態系に及ぼす影響 〜外来植物対策をめぐる課題〜  [B]
W03 地域生態系の保全・再生に関する合意形成とそれを支えるモニタリング技術の開発 [C]
W04 生態から種分化に迫る [D]
W05 ゼロからわかるセイヨウオオマルハナバチ外来種問題 [E]
W06 化学物質の生態影響評価に対して生態学者ができること [F]
W07 病原体と人間社会、生態系を結ぶマルチスケールデータベース構築に向けて [G]
W08 世界自然遺産登録に向けての小笠原研究の現状と展望 [H]
W09 隠すな！ＮＤ（旧総合説）の論理破綻、結集しよう！新今西説を機軸とした“新総合説”へ [I]
W10 第二弾：陸域生態系におけるリン制限仮説と生物の適応 [J]
W11 陸水湿地の環境構造と生物多様性 [K]

◎ 3/20 (TUE)
公募シンポジウム (9:00-11:30) [Room]
S01 ミクロな世界からの新展開： 微生物群集の不均一性が支える物質循環と生物間相互作用 [A]
S02 花粉一粒の直接遺伝解析から見えるもの [B]
S03 人間活動下の生態系ネットワークの崩壊と再生 [C]
S04 大規模開発につける薬（３）――生態学は政治・司法に発言できるか [H]

自由集会 (18:00-20:00)
W12 データ解析で出会う統計的問題 - ベイズ統計学の考えかた・使いかた [A]
W13 群集生態学における適応：進化・可塑性と生物群集の相互関係 [B]
W14 氾濫原の自然再生に向けて [C]
W15 エンドユーザからみたDNA バーコーディング [D]
W16 蜂合わせ 〜蜂好き達の宴〜 [E]
W17 生物の空間分布データから生態学的プロセスを推測する [F]
W18 植物における物質分配の生理生態学 [G]
W19 特定外来生物防除の最前線〜対策と問題点〜 [H]
W20 タケは里山の厄介者か？ 〜竹林の拡大と環境への影響〜 [I]
W21 侵入から四半世紀〜ジャンボタニシは今どうなっているのか [J]

◎ 3/21 (WED)
企画シンポジウム (14:30-17:00) [Room]
T01 Stoichiometric frameworks for evolutionary, population, and community dynamics [A]
公募シンポジウム (14:30-17:00)
S05 長期気候変動と熱帯雨林 〜 熱帯域に於ける森林衰退の背景とその影響 [B]
S06 希少樹種の現状と保全 [C]
S07 海洋性ベントスの生活史：成長と繁殖の生態学 [D]
S08 植物生態学のための山地地形の扱い方． [E]
S09 流域の保全：大規模な止水域が流域環境に与える影響解析 [F]
S10 人間活動がもたらす生態系の変化：常識のウソを探る [H]

◎ 3/22 (THU)
フォーラム (9:00-11:30) [Room]
T02 大規模長期生態学のためのデータベースと情報公開（大規模長期生態学専門委員会） [A]
公募シンポジウム (9:00-11:30)
S11 高山・亜高山帯における生物集団の維持機構―遺伝子流動という視点から― [B]
S12 観測からモデリングへの新たな流れ：地球生態学のデータ利用 [C]
S13 チョウ類とその生息環境の保全をめぐる問題 [D]
S14 京都盆地における環境変遷と人間活動 - 様々な手法によるアプローチ [F]

自由集会 (18:00-20:00)
W22 博物館の生態学３ 〜標本のチカラ〜 [A]
W23 数理生態学はどこからきて、どこへゆくのか [B]
W24 生物の時間空間変動理解のための実証的・理論的アプローチ [C]
W25 淡水産外来無脊椎動物の侵入実態と防除に向けた課題 [D]
W26 樹木のモジュール構造は繁殖にどう関わっているのか [E]
W27 パラサイト・パラダイス〜日本にやってくる侵入寄生生物〜 [F]
W28 群落談話会〜種多様性概念と地域の生物多様性保全 [G]
W29 生態学者よ、街に出よ！Part4 研究活動と企業環境活動の架け橋 [H]
W30 南西諸島における亜熱帯林の動態と生態系保全 [I]
W31 第13 回 Matter Flow and Ecosystems [J]
サテライト集会 種子散布研究会・フェノロジー研究会 (19:30-22:00) 別会場

◎ 3/23 (FRI)
サテライト集会 アグロエコロジー研究会〜 在地・在野の農生態系インタープリター大集合 〜
生態学の学究と農村現場をいかに結ぶか？〜  別会場

• うーむ， べいづなデータ解析自由集会 (W12) と同じ時間帯に またべつのデータ解析っぽい自由集会がかちあってしまっているじゃない …… とはいえ， けっこう「かちあって」る自由集会は多そうに見えなぁ． 自由集会などが多すぎるのであろう． 別の日だけど内容的にかなり似ているのもいろいろあるし．
• とか データ解析サイト ひねくってるうちに， 「R の help() の html 版だけを自分の website に置くと便利じゃないかしらん」 という気がしてきたので， そのあたり取り組んでみる．
• 手もとの ThinkPad にいれてる R package は日々変わるので (R の version だってよく変わるわけだし) …… これは /usr/lib/R 以下の HTML help file を自分の web page 置場 /home/kubo/public_html/ に rsync するしかない．
• rsync は R とは関係ない Unix 系のファイル転送コマンドであり (例によってゐんどーづ版 binary も用意されてるみたいだけど)， たいへんたいへん便利， 私などはもうこれに頼りきって日々の 差分バックアップ・アップロードなどをこなしている．
• 超絶的に便利であると同時に， これは冗談ではなく危険きわまりない兇悪なコマンドであり， 指示をたった一文字まちがえただけで (これも冗談などではなく字義どおり) 「何もかもきれいさっぱりあとかたもなく削除」 できてしまう， という ……
• とゆーことで慎重に小規模な実験を繰りかえしながら R の HTML help の同期をとるコピースクリプトを書く． とりあえず できた． 他にもいろいろと手なおしを …… などとうだうだしてるととっくに午後に突入していることに気づいた． ひと区切りつけて昼飯．
• 1450 自宅発． 午前中は雪ふってたけどいまはほとんどやんでる． 東区の北 9 の Homac (DIY 店) でいろいろと買いもの． あと北大前の TwoTop で予備 HDD を公費買い． タイ製 Segata ST3802110A/72R/2M (Barracuda 7200.9 80GB) 4980 円． 1610 研究室着． 院生密度， 今日も低くはない．
• 予備 HDD のとりつけと設定は …… まあ， 後日にしよう．
• 助教書類かき， やる気がでなくていつまでも終わらず． いやいや， こんなモノに時間かけてちゃダメだな．
• …… とか言いつつ， 作文 2000 字中 700 字で撤退． 2105 研究室発． 2120 帰宅． 体重 71.0kg． 晩飯．
• [今日の運動]
  • うんどう休養日
• [今日の食卓]
  • 朝 (0800): 米麦 0.5 合． マイタケ・ワカメ・ホタテ・豆腐のシチュー．
  • 昼 (1330): うどん． マイタケ・ワカメ・ホタテ・豆腐のシチュー．
  • 晩 (2200): 米麦 0.7 合． 卵丼．

2006 年 12 月 24 日 (日)

• 0750 起床． 静かな朝だな …… と思ったら昨晩のうちにかなり積雪してた． つもった新雪による消音か． いっきに札幌の冬らしく． 朝飯． コーヒー． 洗濯． 怠業．
• 1100 自宅発． 豪雪． JR 札幌駅内うろうろしてから 1140 研究室着．

[いきなり豪雪]

目測で 40-50cm ぐらいいっきに積もったかんぢ． 少ない雪に困ってた札幌市周辺スキー場は大喜びだろう． おそらく朝の短時間につもったのか， 除雪がぜんぜん追いついてない． 札幌駅周辺はハンパな融雪システムによって泥濘化． 北大につくころにはなぜか晴． 地環研にいたる経路は雪にうもれてたけど．

• また web site 上の R help html ファイルの CSS ファイルその他いろいろをひねくってみる． ユーザー数最大の Microsoft Internet Explorer 6.0 だけ が他のブラウザーとは異なる CSS 解釈 (W3C 勧告無視の MS オレ流手前勝手な解釈) するので毎度のことながらつくづくイヤになる． 今回は対策せずに放置． 昼飯．
• やはりちょっと気になるので調べてみたら …… なンと MSIE6 は DOCTYPE 宣言で 挙動が変化するのか! …… しかもその状態遷移は奇怪そのもの． こんな阿呆阿呆な MS 仕様は無視したいところだが， やはりなんだかとても気になってしまったので …… たとえば「ぎょーむ日誌」なんかだと

  <!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN">
  

  だったのを

  <!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN"
   "http://www.w3.org/TR/xhtml1/DTD/xhtml1-transitional.dtd">
  

  と URL をつけてみる …… MSIE6 における表示が変化してしまった． MSIE6 の動作が MS 僭称する「標準準拠モード」に移行したため …… 何が「標準準拠」だ， ちくしょうめ． 敗北感あるなぁ． いやはや． やれやれ．
• どーりで Strict なペイジを VikiWiki に生成させてる 生態学のデータ解析 サイトはレイアウトのくずれが MSIE6 でもほとんどなかったわけだ ……
• 投げやりなきもちがつのってきたので， 「R の help ファイル」 公開用に ~/public_html/Rdoc/ 以下においてる全 *.html ファイル (現在の総数 4913 個) 中の 「Transitional なる文字列をことごとく Strict に置き換えよ」 という乱暴な命令:

  find ~/public_html/Rdoc/ -name '*.html' | while read html;  \
  	do perl -i -p0777e "s/Transitional/Strict/o" "$html"; \
  done
  

  を発してしまいましたよ (処理時間およそ 20 秒)． 敗北感もますますつのる．
• 敗北感かみしめつつ， お茶部屋で院生たちとお茶 & ケイキ．
• R の help HTML ファイルは， そもそも package インストール時に Rd ファイル (R のドキュメント形式) から /usr/lib/R/share/perl/build-help.pl なる Perl スクリプトによって生成されたものであり， DOCTYPE 宣言はその子分 Perl モジュール内できめうちにされている …… と判明した． 改造はやめといたほうがよさそう．
• いやいややってた 助教書類 かき (次の水曜日〆切)， 「いつまでもこんなことにだらだらと時間とられていてはダメだ」 と， きあいで 2000 字作文を終わらせる． しかし時刻はすでに 1755．
• ちょっと作文その他のみなおし． まあ， 明日もういっぺんみなおしてさっさと提出してしまおう …… ということで本日は終了． 1850 研究室発． 積雪でおおわれた札幌市内ふらふら散歩． 晩飯． 帰宅．
• [今日の運動]
  • うんどう休養日
• [今日の食卓]
  • 朝 (0830): 米麦 0.5 合． 海藻スープ．
  • 昼 (1420): 研究室お茶部屋． 食パン．
  • 晩 (1940): 大通公園横のスープカレー屋， Spice Labo．

2006 年 12 月 25 日 (月)

• 0750 起床． 0810 自宅発． 晴． 0825 研究室着． 朝飯． コーヒー．
• とある共同研究の文書よみとそのメイルかき．
• そしてそのあと正午まで， お茶部屋で宮田さんに先日の苫小牧 13 樹種のシュート内重量分配 & SLA の階層ベイズモデルに関する説明．
• つぎに庄山さんと知床の樹高成長モデリングに関する相談 …… というか， 二週間ほど前 に

  今回は 「モデル・計算がちょっとめんどうになったら Bayes + MCMC 計算」 原則を適用して WinBUGS にめんどうを押しつける予定．

  とか書いてたあの一件で， 私は見とおしの甘いてきとーなる自分自身の発言の責任 とらねばならぬ状況に追いつめれてるわけだ． 近ごろ， こういうパターンが多いような気がする． だれか助けてー
• 昼飯． お茶部屋ではコンヴィニエンスストアのクリスマスケイキの大人買い (半額になったので丸一個購入) & 子供食い (一人でそれ全部を食べること) とよばれている 珍しい現象が観測された．
• 院生たちの年末雑談傍聴．
  「部屋の大そうじがたいへんでさー」
  「そうそう， 帰省で飛行機のる前にきちんとそうじしなきゃいけないもんね」
  「えっ? どうして?」
  「だって飛行機とか落ちて死んだときに， 遺品の回収きたヒトたちにちらかった部屋みせるわけにはいかないぢゃん． って， 親にそう言われた」
  「……」
• またあちこちメイル書き & データ解析雑談．
• 昨日ひととーり書いた助教書類のみなおし． 「今後の研究・教育についての計画と抱負 (2000 字程度)」 ってところに 「環境科学における統計学的手法とその教育の重要性」 と題して何だか地環研生物系の統計学的手法教育の 貧弱さをねちねちと指摘した作文をものしてしまっている． こんなので大丈夫なのだろうか? …… よく考えぬまま送信． ともあれ一件落着．
• 次の〆切つきぎょーむは 1/9 (火) 夕方の 生態学会松山大会 講演要旨の登録 か．
• 年末年始の楽しい 節句働き 予定一覧 (順不同):
  • 知床樹高成長 Bayes モデル
  • 石狩浜ハマエンドウデータをゑくせるファイルから救出
  • サロベツやちパワーモデルの計算やりなおし
  • 母子里林冠 MCMC 計算の投稿論文原稿の修正あれこれ
  • 大雪山系地衣類繁殖モデルの作文
  • (松山大会一般講演の) 屋久島葉齢分布 Bayes モデルの計算やりなおし
  • A 棟 8F ほか越年隊の皆さんのデータ解析こんさる各種
  • ……
  うーむ ……
• 「放置するとコワい順」を慎重に比較検討 …… いくつかは比較不可能なほどにコワいんだけど， とりあえずまずは知床問題から着手． 図を何枚か作図したところでばててきたので撤退． 1810 研究室発． 1825 へろへろと帰宅． 体重 70.4kg． 晩飯の準備．
• 厳冬期用防寒具 (Armée Française M-64 Parka, カラー： OD, サイズ： 92L = 身幅 71cm 着丈 97cm 裄丈 92cm) を宅配便で受けとる． もちろん今まで使ってたやつも現存してるんだけど， すでに咋シーズンの段階で 「なンでそんなボロボロになってるものをいつまで着てるんですかっ!」 「そんなの着てこのへんをうろうろしないでください」 といった A 棟 8F における 政治的圧力が抵抗しがたいほどに強まってしまっていたので ……
• 晩飯． 晩飯くったらいきなり眠くなり， あやうく椅子の上で昏睡しそうになった． 買いものに出て (と言ってもとなりの北 12 条生協だけど) ふらふらしてるうちに眠気がさめる．
• [今日の運動]
  • またうんどう怠業日
• [今日の食卓]
  • 朝 (0830): 研究室お茶部屋． 食パン．
  • 昼 (1300): 研究室お茶部屋． 食パン．
  • 晩 (1950): 米麦 0.7 合． マイタケ卵丼．

2006 年 12 月 26 日 (火)

• 0600 起床． 朝飯． コーヒー． 0800 自宅発． 曇． 0815 研究室着．
• 知床樹高成長のベイズモデリング …… これって無情報事前分布ばかりで階層ベイズになってない． というのも樹種を Abies とそれ以外， として計算してるんで．
• いや， そんなことない． 観測誤差の事後分布は超事前 (hyper prior) 分布 から生成させているから階層ベイズモデルだ．
• 1030 より 研究室セミナー， 本日は西村さんで サロベツ湿原調査地 で施肥実験 (窒素肥料投入) やりました， というもの． どうもハナシとしてはサロベツなんかでは泥炭採掘とかで湿地を堀りかえして 植物がいないときが地下水中の窒素濃度が高く， 植物がはいってくるとこれが下がるらしい． 「極相」ということになってるミズゴケ群落では窒素が少ない． で， 窒素肥料を追加投入してやると， なぜかしらミズゴケなんかでは群落高とかが低くなったりする (原因不明， 「肥料まけ」説もあり)． 外国での研究例なんかだと， 大気汚染だか地球環境変化だかによって雨中の窒素濃度が高まり ミズゴケ群落がこれによって影響をうけてしまって …… というハナシになるらしいんだけど， サロベツではこれはいったいどういうハナシになるのかは まだ考え中とのこと．
• 知床計算ひととーりできた． Abies (トドマツ) は収束いいけど， 広葉樹 (個体数少ないのでまとめたもの) はダメだ …… とりあえず発注者の庄山さんに「ここまでできました」 説明． じつはこれは自分のアタマの整理に有用で， 「広葉樹のパラメーターはまとめたままだけど， 定数部分に (広葉樹間の) 種差の random effects いれてみるか」 というアイデアを思いつく．
• またいろいろ計算． 昼飯． 帰省ラッシュ・高額航空券を回避した院生たちの移動により， 院生密度がぢりぢりと低下中．
• 年末年始でおもいだしたんだが …… 来月 20 (土) & 21 (日) 日のセンター試験監督とかってやらなくてもよいのかしらん? 過去のぎょーむ日誌によると「召集令状」は毎年 11 月末 - 12 月始めにでてたらしいんだけど， 今回は私のところには連絡ナシ． 4 回連続でやったら 1 回休みにしてもらえるわけ?
• 知床 Abies + 広葉樹 (合計 19 樹種) の樹高成長モデルに 「(広葉樹間) 種差 random effects」 を導入して， これを例によって R2WinBUGS で MCMC 計算 (Gibbs sampling)， という作業つづく． 広葉樹間の種差なし DIC = 3269 で 種差あり DIC = 3265 なんで， びみょーだけど種差ありでやりますか (トドマツとそれ以外の間の種差は最初から考慮している)．
• 結果はこんなかんぢで． これはトドマツ (Abies) 稚樹の高さ成長の 4 パラメーター， 定数部分 [1]， log 樹高依存性 [2]， シカよけ柵 (exclosure) 処理効果 [3]， カラマツ林冠木まびき (thinning) 処理効果 [4]． 一般化線形モデル的な階層ベイズモデルでよく説明できていそうに見える．

• (4-5 年前の?) カラマツ林冠木まびきでトドマツの成長が悪くなる理由は， 広葉樹稚樹のほうが応答はやくトドマツ稚樹を被陰するから， と庄山さんの説明． 結果をいろいろと説明してから， ファイルあれこれを庄山さんにひきわたす． これにて年末節句働きのひとつ， 一件落着 …… か?
• またど忘れしないうちにここにメモしとくか …… 上の樹高成長モデルは 平均 exp(線形予測子) のガンマ分布 (分散 ∝ 平均， 試してみたけど 分散 ∝ 平均² はダメ) が「真の」高さ成長量で， それに観測誤差が正規分布で追加されるという 「浦口さんカエデ成長モデル」 のベイズモデル版であります， と．
• 4 年前の今ごろ (このころは R2WinBUGS はおろか R もマトモに使えなかったので Perl で計算 + gnuplot で作図 という時代) であった浦口さん計算のときはガンマ分布 * 正規分布たたみこみの積分をふくむ 尤度方程式の最大化 (つまり最尤推定) で観測誤差の大きさも点推定していた． 今回は (R2WinBUGS から呼びだす) WinBUGS でもって事後分布を推定 (「モデリングや計算がめんどうになったらベイズ」ぽりしー)， 最大事後確率 (maximum a posteriori; MAP) 推定値だと， 稚樹たちの樹高成長量の観測誤差の標準偏差は 8-9cm ぐらい．
• 次． われながらヒドいことに 先月 27 日 いらい (MCMC 計算の収束がうまくいかなくて) 放置しているサロベツやちパワー計算の再開を画策． まずは小山さんのところにいって平身低頭 & 現状の確認と今後の相談．
• 事務からまたおかしなメイルが …… なにー? また「助教」書類だせって? 先月の命令 で作文を昨日送信したばかりなのに， また同じような 新しいやつ が …… これって何がちがうのかしらん? わけわからんなぁ …… まあ書類だせってのは 1/15-1/31 なんでしばらく放置するか．
• うわー， 香川大の小林さんから「もうひとつの」アカマツ原稿が ……
• サロベツ計算をススめなくては …… と思いつつも急に空腹感が自覚されてしまって …… 撤退． 1840 研究室発． 1855 帰宅． 体重 70.4kg． 晩飯の準備． 晩飯．
• うーむ， 夜は仕事がススまない ……
• [今日の運動]
  • 腹筋運動 30 × 3 回． 腕立ふせ 5 × 3 回．
• [今日の食卓]
  • 朝 (0630): 米麦 0.5 合． ネギ卵炒飯． 海藻スープ．
  • 昼 (1310): 研究室お茶部屋． 食パン． リンゴ．
  • 晩 (2030): 米麦 0.7 合． キャベツ・ニンジン・ブナシメジ・コンニャクのシチュー．

2006 年 12 月 27 日 (水)

• 0855 起床 …… えーと， 夜明け数時間前になぜかめざめて， 「うー，眠れん」 状態がつづいて， 夜明け後にまたちょっと寝てこういう時刻に起床しました． 年末ぼけかな． 0910 自宅発． 南方からやってきたヘンな低気圧のせいで異様に気温たかく (本日の最低気温は氷点下ではない!)， 雨． 0925 研究室着．
• お茶部屋で朝飯 & コーヒー …… を準備してたら庄山さんにすばやく捕捉され， 昨日の知床樹高成長の階層ベイズモデル MCMC 計算 に関する質疑応答あれこれ， に一時間ほど．
• 宮田さんから苫小牧 13 樹種データに訂正ありました， とかで新データファイルが送られてくる． これって「計算やりなおせ」ご命令だろうなあ …… ということで， ssh で Dell Linux 機に MCMC 計算やりなおしを命じる．
• 一連の行列計算でケリをつける (やや複雑なる) やちパワー計算用の小山さん R プログラムの解読・試験運転にとりくむ．
• とりあえず私の環境でも今までの計算は再現できるようになってきた． 昼飯． 院生密度は順調に低下中．
• すでに内地方面への機中のヒトとなっているであろう宮田さんからの 「宿題」 MCMC 計算， 6500 秒を費して終了 (やはりデータが変わってもこれぐらいの計算回数が必要だった)． うーむ， 前回 とびみょーに違うな． データ数の変化は +1． 前回 「どう考えてもおかしい」 として外した一個体の観測値が修正されたため．

• まあ， よくわからないのでこれは放置．
• サロベツ湿原上の植物分布がやちパワーに依存してるかどうか問題， 単純に計算すると「やちパワー」が最大になってしまうやち坊主の上には あまり植物がいない (やち坊主を形成しているスゲのシュートは高密度に存在するけど)． このあたり， どうモデルにくみこむか …… 計算の単位になってる 10×10cm 区画が「やちの上にいるかどうか」を判定する関数を追加， これを説明変数としてみる， というのはどうだろうか． とりあえず 1 line 内 1000 区画について計算し， 6 line ぶんの図を作ってみてこの計算結果をチェックする．

> apply(c("Non-Top" = 0, "Top" = 1),
     function(tt) summary(as.factor(d.grid[d.grid$Ttop == tt, "Freq"])))
$`Non-Top`
     0      1      2      3      4      5      6      7      8     10     12     15 
122143   1837    283     69     29      6      6      3      2      1      2      1 
    16     19 
     1      1 
$Top
   0    1    2    3    4 
5528   63   17    5    3 

• これを例の 「WinBUGS でも理解できる ZIP モデル すなわち隠れ変数版 ZIP」 に組みこもうとする試行錯誤．
• ひととーりできたので， ばぐ取りにかかろうとしたら， ちょっとした server 頼まれごとが． そのご褒美に珍しい氷河写真ファイルをいただいてしまった． わーい

[Laguna San Rafael 国立公園 (拡大)]

チリのやや南方 Aysen 地方のフィヨルドにつきだした San Rafael 氷河の末端 …… 今月はじめに撮影したものだそうです．

• サロベツ改良モデル， Solidago (アキノキリンソウ) を実験台につかって計算プログラムのばぐとり．
• 北村君と阿寒アカエゾマツ成長解析のこんさる少々． lmer() でケリをつけることになりそう ……
• サロベツは全優占種に関する種ごとのベイズモデル MCMC 計算 (そう， このサロベツ問題は種差 random effects 法を適用できていない!) を命じておいて撤退． 1955 研究室発． あいかわらず， 雨． 2010 帰宅． 体重 70.0kg． 晩飯．
• やはり夜は何もススまない．
• [今日の運動]
  • なのにうんどう休養日になりがちな年末
• [今日の食卓]
  • 朝 (0950): 研究室お茶部屋． 食パン．
  • 昼 (1310): 研究室お茶部屋． 食パン．
  • 晩 (2100): 米麦 0.7 合． キャベツ・ニンジン・ブナシメジ・コンニャクのシチュー．

2006 年 12 月 28 日 (木)

• 0530 起床． なぜかしら今朝もふつーよりかなり早い時刻に目がさめた． コーヒー． 歳末 R たすけあい運動実施中なので， 朝からとあるデータ変換プログラミングを楽しむ．

• 朝飯． 0900 自宅発． 曇． 昨日の雨でかなり雪がとけてしまった． 0915 研究室着．
• サロベツ湿原調査地 やちパワー計算， Moliniopsis まで終了して Lobelia で Trap 窓だして計算がコケていた ……

> v.spc.count[v.spc.count > 100]
Rhynchospora  Hypochaeris      Drosera     Solidago  Moliniopsis      Lobelia   Eriophorum        Carex 
         886          518          372          350          306          173          123          112 

• そしてあいかわらず MCMC 計算の収束はよろしくなかったり． さすがに ひと月前 よりはマシだけど． ああ， でもこれってモデル改善というより， (私があまり好きではない) 「説明変数中央化 (X <- d$X - mean(d$X)) MCMC 計算こざいくワザ」 の影響かもしれん．

• あー， 改善策がおもいつかない ……
• よくわからないので超事前分布， あのよくわからぬ Wishar 分布に与える対角行列を diag(0.05, n.r) から diag(0.01, n.r) と変えてみる． また Solidago で試験運転 …… 収束の良くなさも含めて， ほとんど何もかわらん．
• このあたりの挙動忘れてしまったので， diag(1.0, n.r) で再計算 …… やはりあまり変わらん．
• 超事前分布として使ってる Wishart 分布のパラメーターを diag(0.05, n.r) にもどして …… 次の一手， 「『やち (坊主) 上』効果は Poisson に限定」 モデルはどうだろう? 「植物がいない」ことを説明する部分を複雑にしすぎるとダメだろう， という気がしたんで． これはめちゃくちゃなコトを言ってるみたいだけど， 必ずしもそうではない． ZIP 的なモデルは「いる・いない」 + 「いるとすればどれぐらい?」 からなりたっているんだけど， 観察される現象は
  1. やち嫌いの植物: やち上だろうがやち周辺だろうが， やちパワーの高いところはともかく避ける
  2. やち好きの植物: やちパワーの高いやち周辺は 存在確率 & 密度が高い， ただしやち上はかなり密度が低い
  3. 極端なやち好きの植物: やち周辺・やち上どちらでも 存在確率 & 密度が高い
  となっているので， 「やち上かどうか」がモデルの中で使われるのは 2. の場合だけなのである． このときにやち上の区画で「(仮に) 存在確率は高いが，密度は低い」 となっていれば現象を説明するのに不都合はない． あ， 「存在確率だけが低い」 としてもよかったか． うーむ．
• ともあれすでに計算はスタートしてしまっているのでそのまま続けることに． 30 分ほどかかる計算まち時間のあいだに昼飯調達にでる …… といっても， 今日はめんどうになったので地環研ロビーの北大生協の出張弁当販売で．
• 単純化したモデルの計算終了． さすがに収束はマシになった． もうちょい単純化できないかしらん?
• 存在確率および密度の 「切片」に該当する部分を無条件事前分布から生成させてきたんだけど， これをやめてみたらどうだろう? つまり平均ゼロの無条件事前分布にすべてをまかせてしまう， と …… お茶部屋で昼飯くってるうちに計算が終了したんだけど， うーむ， これはいまいちだな．
• ということでもとにもどして， 今度は存在確率が「やち上かどうか」に影響されるようにしてみる …… まあ， こちらのほうが「説明」はやりやすいだろう．
• 次の単純化． これまで stage (植物の発育段階) によって「やちパワーのききかた」が違う， としてたんだけど …… この観測データは (意外と情報量が少ないため) 複雑な統計モデルをもちいて情報を抽出できない， ということが次第にわかってきたので， 「どの stage でもやちパワーの影響は同じ」 としてみる (というのもこのあたりのパラメーターの収束は劣悪なんで) …… Solidago のみの試験運転の結果なんだけど， 予想どーり収束は格段によくなり deviance は改善されている．
• といった状況を小山さんに説明して， このサロベツやちボスの了承を得て主要 8 植物種に関する計算を始める． 予想される計算時間はおよそ 5 時間 (途中でコケなければ)． 現在時刻は 1410．
• サロベツつかれを感じて (歳末効果ですっかり静かになった) お茶部屋で静かかつくたびれお茶部屋雑談．
• 「年末年始の楽しい節句ばたらき一覧」 とやらをみなおしてみる．
  • 知床樹高成長 Bayes モデル
  • ひとまず終了?
  • 石狩浜ハマエンドウデータをゑくせるファイルから救出
  • まー， 稲葉さんはいま九州 → 中部 → 関東周遊の大旅行中だし ……
  • サロベツやちパワーモデルの計算やりなおし
  • 現在計算中 …… よーやく終わりが見えてきた?
  • 母子里林冠 MCMC 計算の投稿論文原稿の修正あれこれ
  • まだ何も……
  • 大雪山系地衣類繁殖モデルの作文
  • まだ何も……
  • (松山大会一般講演の) 屋久島葉齢分布 Bayes モデルの計算やりなおし
  • まだ何も……
  • A 棟 8F ほか越年隊の皆さんのデータ解析こんさる各種
  • 急患随時うけつけ・処置中
  ということで， いそぎの用事もないし， 小林さんから先日おくっていただいた 例の「もうひとつの」アカマツ論文原稿にとりくんでみるか．
• えーい， なかなかアカマツ頭にきりかわらん．
• 1859 サロベツ計算終了． 1200 MCMC step 8 植物種におよそ 4 時間 40 分， か． そして結果は今まで得られた中で (さすがに) 最良であるけれど， まだまだ収束よろしくなく， そして不思議な事後分布が含まれている (まちがいではないのかもしれないけど)．

• もっと長い (2000 step) 計算を命じておいて撤退． 1930 研究室発． 晩飯． 2055 帰宅．
• [今日の運動]
  • またうんどう休養日
• [今日の食卓]
  • 朝 (0710): 米麦 0.7 合． タマネギ・卵の炒めもの．
  • 昼 (1250): 研究室お茶部屋．
  • 晩 (2010): 北 11 西 1 の「でみ蔵亭」， 焼きハヤシ (ハヤシライスのグラタン)．

2006 年 12 月 29 日 (金)

• 0950 起床． うう， いきなり年末年始休業もーど， か …… 朝飯． コーヒー． 洗濯． 怠業． 1130 自宅発． 曇． 今日も気温高い． 現在， +1°C ぐらいか? 1150 研究室着．
• 昨晩帰るまえに計算させていた 2000 MCMC step のサロベツやちパワー ZIP ベイズモデル， またまた Lobelia (サワギキョウ) で WinBUGS がコケてしまって， Moliniopsis (ヌマガヤ) まで (優占 8 種のうち第五種まで) で計算が終了していた． とりあえず， 重要なパラメーターの事後分布の作図．
• うーむ， これは時間のかかる計算だったというべきなのか …… とりあえず上の図のパラメーターはだいたい収束しているようにみえる．
• 小山さんとやちパワー地図などと照合しつつ 上の結果を検討してみたんだけど， だいたいまあこんなもんでしょう， と． ミカヅキグサはやや意外なかんぢだが， 他のデータとはたしかに整合しているとわかった．
• 昼飯． 院生密度はさらに下降． ところで， 今日はすでに大学としては年末年始休みなんだけど， 地環研の全館暖房スチームヒーターは作動している． これは冬でも土日などは止められてしまうのに， 何で今日はそうなってないのか? 露崎さんによれば凍結防止だそーで． ふーむ． 三日以上とめるとアブない， とか?
• Lobelia 問題， 43 分目でコケる， という現象が再現された． 2000 MCMC step 計算したあとの DIC 計算とかそのへんかしらん? よくわからんのでパラメーター初期値を変えて再計算を命じる． とりあえず優占種グループ上位 5 種の計算はすでに終了しているので， (重要度のやや低そうな) 優占種グループの下位 3 種はまあそんなにあわてなくてもよいわけで．
• 矢澤さんに Pasoh 熱帯林 50 樹種の太さ-高さ関係の最新データをみせてもらう． これって意外と allometric な関係が成立しているようにみえたりする …… つまり「アタマうち」とかが明瞭ではなさそう． ふーむ． そして宮田さん苫小牧 13 樹種シュートモデリングが それなりに成功した現在においても D-H 関係はどうモデリングしたものか， じつのところよくわからない． 「材積」を (制約条件的な) latent variable にして， といった妄想もあるンだけど ……
• サロベツ Lobelia 計算， 初期値かえたら今度は開始 10 分後にコケた． うーむ． まったく挙動不明なのでまた初期値をかえて再計算発令 …… 今度は 13 分後にコケた
• なんかサワギキョウのデータってヘンなところがあるにちがいない， と思って調べてみると …… あった． 全 6 調査 Line 上に必ず出現するにもかかわらず，

  > tapply(dL$Freq, dL$line, sum)
    A   B   D   E   G   H 
    6   1  15   3 134  14 
  

  全 173 個体のうち 134 個体が調査 Line G に集中してるんだよなぁ …… じつに極端なる集中分布だ． たぶんこのあたりに原因ありそう．

• ちなみに優占 8 植物種の Line 間ちらばりはかくのごとし． ブタナ (Hypochaeris) みたいに A-E Line に存在しない場合には， Line ごとデータからぬかれているわけだが ……

  > sapply(v.target.spc, function(s) { d <- dall[dall$spc == s,]; tapply(d$Freq, d$line, sum) })
  $Rhynchospora
    A   B   D   E   G   H 
  183  20 298 268 104  13 
  
  $Hypochaeris
    G   H 
   74 444 
  
  $Drosera
    A   B   D   E   G   H 
    3   9 160  33 123  44 
  
  $Solidago
    A   B   D   E   G   H 
    3   6  63   4 230  44 
  
  $Moliniopsis
   A  B  D  E  G  H 
  97 81 31 39 56  2 
  
  $Lobelia
    A   B   D   E   G   H 
    6   1  15   3 134  14 
  
  $Eriophorum
   A  B  D  E  G  H 
  40  9  9 19 14 32 
  
  $Carex
   A  B  D  E  G  H 
  34 15 10  9 27 17 
  
  

  まあ， 「1 Line だけ集中」 という Lobelia は変わりものなのかもしれない． しかも， 特定の stage に偏っているし，

  > tapply(dL$Freq, dL$stage, sum)
    f   j   s 
    8 153  12 
  

  ということで， とりあえず 個体数 5 未満の Line は計算に含めないよう改造して， また再計算発令 …… また 13 分でコケたので， 条件をキビしく個体数 10 未満は捨て， としてみる …… これでもダメか …… ならば Line G だけ， ならば? …… うーむ， これもダメ．
• 歳末急患的 R こんさる． vector の中である条件をみたす値だけを別の vector の値で置き換えたい， すなわち「条件つき vector コピー」はどうやったらよいのか? R はこれを実現する独特の文法をもつ:

  > dbh.new <- c(10, NA, 21, 16, NA, 13)
  > dbh.old <- c(9,   8, 20, 15, 12, 11)
  > (replacement <- is.na(dbh.new)) # どこが欠側値か
  [1] FALSE  TRUE FALSE FALSE  TRUE FALSE
  
  > dbh.new # 欠側値を含む vector
  [1] 10 NA 21 16 NA 13
  
  > dbh.new[replacement] <- dbh.old[replacement] # 欠側値のみ置き換え (条件つきコピー!)
  > dbh.new # 欠側値を置き換えられた vector
  [1] 10  8 21 16 12 13
  

  …… と書けばよい． 条件式 target は TRUE, FALSE の vector (つまり logical 型) であれば何でもよく， is.na() とかじゃなくてたとえば

  > replacement <- dbh.new < 10.0
  

  などでも良い．
• サロベツ計算は 「データぬきわざ」をやめて「MCMC 計算の初期値工夫」 の方向で． なんとなくうまく動いてくれそうな初期値セットらしきものを さぐりあてたので， 計算を命じておいて本日は終了． 1940 研究室発． 北 9 東 5 の Homac 光星店で後述する「とあるもの」を購入． 2030 帰宅． 体重 70.4kg．
• さて， 私はここで自分が長年にわたり 「ダンボールに寝る生活」 を続けていたことを告白せねばならぬ． 2001 年 3 月 11 日に現在の「高床式」ベッドを入手したんだが (ぎょーむ日誌 目次 でいうところの「愚行の柏」事件， 当時は東京都大田区大森東に住んでた) …… ほんらいの使用にはマットレスが必要とされるこのベッド， ダンボールをしいてそのうえにふとんをのべていたりしました． はい．
• しかしベッドの上のダンボールたちも林床に成層するリターのごとく， しだいしだいに圧縮されるみたいで (しかも関西空港と同じく「不等沈下」する) …… 近ごろはいささか寝心地よろしくなく (腰のあたりが痛い)， そのせいか「朝になってないのにめざめる」ことが頻発， そして新しく追加すべき手もちのダンボール箱も尽きてしまったので， 今回のマットレス購入にいたった次第である． お茶部屋の新聞のおりこみ広告にちょうど Homac の安売り情報がまぎれこんでいて， ですね．
• ということで， さっそく交換作業にとりくむ．

[ダンボール「マットレス」]

しきふとんを除去した状態． このようにダンボールたちが成層しております． このベッドは「脚」が長く 1m ぐらいある …… ので， 下に荷物をおけたりする．

[ダンボール除去]

ダンボールを取りのぞくと骨格とグリルが …… のはずなんだけど， 腰にあたる部分に「補強板」が入っているのを発見して驚愕する． そうだった …… ダンボールだけでは「強度」が足りなかったんで (腰のあたりだけが沈下していく)， こういう応急処置したんだっけ． どうりで寝心地が悪いわけだ．

[マットレス展開]

我が家から 1000m ばかり離れている Homac から人力輸送してきた (頭上に乗せて運ぶ方式がラクでした) マットレスを伸張． もちろん「補強板」は除去ずみ． 脱ダンボール寝 はじつに快適で， 少なくともこの晩は安眠できました …… なんだか「帰ってきたんだ」感ありました．

• 晩飯の準備． 晩飯．
• [今日の運動]
  • またまたうんどう怠業 ……
• [今日の食卓]
  • 朝 (1020): 米麦 0.6 合． ネギ卵炒飯．
  • 昼 (1340): 研究室お茶部屋． 食パン．
  • 晩 (2230): 米麦 0.7 合． 卵焼き． キャベツ・ネギ・豆腐・煮干の味噌汁．

2006 年 12 月 30 日 (土)

• 0940 起床． 脱ダンボール寝の快適さについつい長時間寝てしまいました． 朝飯． コーヒー． 外は雪． また積雪している．
• 怠業 …… しつつ樹木の形状に関するベイズ的な統計モデルの検討． これまでのそういう研究 (屋久島葉寿命・ アカマツ軸成長・ 苫小牧 13 樹種シュート内資源分配など) から観測値の背後にある「真の値」(仮称) の 「可動範囲」 をさらに制約してやる隠れ変数 (latent variable) をうまく設定設計してやることが， こういった種類のモデリングにおいて重要らしい， とわかってきた (これは私独自の考えとかではなく， J.S. Clark 御大とかが Bayesian 総説なんかで ややわかりにくくあれこれと述べていること)．
• で， 考えるのがなかなかめんどうなのでイヤな樹木の太さ-高さ 問題の場合だと， 「見えない状態」による制約として使えそうなものは何か? いろいろあるんだろうけど …… たとえば材積とか樹齢とか， かしらん? たとえば樹齢 (連続値) の場合なんかだと， 事前分布を平均 100 とかの指数分布にして， …… といった工夫が必要だろうな． さらに何かうまい超事前分布とかを設定できるだろうか?
• 考えまとまらぬまま 1145 自宅発． 雪． 1200 研究室着． お， 今日も全館暖房．
• サロベツ計算， WinBUGS に巧妙な初期値を与えることで難所 Lobelia とその次の優占種 Eriophorum の計算はそれぞれ一時間ほどで終了していた． しかし優占 8 種の最後， Carex はコケていたんだよね． 理由はおそらく Lobelia に最適化した極端な初期値 (きわめて低密度) がまずかったんだろね．
• とりあえず計算終了した Lobelia (サワギキョウ) などのパラメーター事後分布を図にしてみてみる． これも以前の計算とだいたい同じで， 存在確率を決める logistic 部は「やち (坊主) 好き」 なんだけど， 局所密度をきめる Poisson 部は「やちと距離をおきたい」， 背反したものになっている．
• これは推定計算まちがいのたぐいなのかというと， どうもそうではないらしく …… 下は昨日も示した Lobelia のやちパワー地図なんだけど， とくに (高密度な) Line G においてやちまわりが好きなのにやちとは少し距離をおいてみたい， という状況が図示されているように思える．

• さて， 最後に計算しくじった Carex の再計算だけど ……

  > sapply(v.target.spc, function(s) { d <- dall[dall$spc == s,]; tapply(d$Freq, d$line, sum) })
  (略)
  $Lobelia
    A   B   D   E   G   H 
    6   1  15   3 134  14 
  $Eriophorum
   A  B  D  E  G  H 
  40  9  9 19 14 32 
  $Carex
   A  B  D  E  G  H 
  34 15 10  9 27 17 
  

  これまた昨日も示したこの Line ごとの分布状況からみて， Carex の Line ごとの ``random effects'' の初期値はゼロ (つまりありがちな初期値) にもどしてやれば問題なさそうだ． 再計算命じる． 一時間ぐらいかかるだろう．
• 昼飯． 昨日が修論第一稿の「露崎さん〆切」だったので， 本日の院生密度はますます低く ……
• Carex 計算ぶじに終了． これで今年のサロベツ計算は終了だろうか? とりあえずパラメーター事後分布図．
• そしてまた A802 の CANON LBP-2810 の Magenta トナーぎれ …… というか残量 50% なのにトナーがでなくなりつつある． これまた不良な再生トナーカートリッジをつかまされてしまったようで． 無料交換してくれるからいいんだけど， 再生屋はこれで商売が成立しているのか?
• なんだか今日はいつもにもましてのんびりとした作業進捗． 窓の外は雪．
• そして「とびこみ的」というか北村君の 阿寒アカエゾマツ年輪データから成長モデルを構築するのに数時間を …… 時系列データつまり一個体から何度もデータとってることになるんで random effects の考慮は必要不可欠， それをなんとか R の lmer() でケリをつけよう， という試み． さらにめんどうなのは「短期的」個体差 (樹木個体内のまいぶーむ?) も同時にモデル化したいんだけど …… 本来なら latent variable を使って Bayesian にして MCMC 計算して， という流れにもちこみたいところだけど， 時間的制約からしてそれは無理． まあ， やや不本意なモデリングになってしまい， しかもモデル選択の結果， やや複雑な構成のモデルが選ばれてしまった． うーむ．
• とりあえず修論ぐらいには使える解析結果が得られたので， 本日は撤退． 2140 研究室発． 雪． 2205 帰宅． 晩飯．
• [今日の運動]
  • 年末うんどう怠業もーど
• [今日の食卓]
  • 朝 (1010): 米麦 0.5 合． キャベツ・ネギ・豆腐・煮干の味噌汁．
  • 昼 (1330): 研究室お茶部屋． 食パン．
  • 晩 (2250): うどん． キャベツ・ネギ・豆腐・煮干の味噌汁．

2006 年 12 月 31 日 (日)

• 0740 起床． 朝飯． コーヒー． 怠業． まあ， ちょっとは歳末らしく， この汚らしい部屋を憲法が保障するところの 「最低限の文化的な生活」 水準まで向上させられないか， とちょっとばかり片づけ & 掃除してみる努力 …… ちょっとばかりの努力でどうこうできる問題ではないな， というサトりの境地に到達できたので， 1120 自宅発． 晴． まだねばっている院生がいるかも， ということで北 12 生協で差し入れ用ミカンをハコ買い． 1140 研究室着． お， 今日もスチームヒーターが作動している．
• web server 機に予備 HDD (80GB) をとりつける作業． 装着そのものは 15 分ほどで終了した． 次に主 HDD が壊れたときに予備 HDD に切り替えるため 参考になりそうなことあれこれ．
  • 新しく買った予備 HDD はケイブル位置で Master/Slave を指定するタイプ
  • しかし Primary の Slave に指定するのは位置関係上 めんどうだったので (HDD 入れ換えもめんどう)， ふだん使わない CDD の IDE ケイブルをぬいて 予備 HDD につけた: つまり予備 HDD は Secondary の Master であり /dev/hdc である
  • えーと， /dev/hda の現状はかくのごとしなので ……

    $ df -h
    ファイルシステム    サイズ 使用中 空き 使用% マウント場所
    /dev/hda1              18G  785M   16G   5% /
    /dev/hda3              19G  6.9G   10G  39% /home
    none                   89M     0   89M   0% /dev/shm
    

    /dev/hdc の partition table 設定

    Disk /dev/hdc: 80.0 GB, 80026361856 bytes
    255 heads, 63 sectors/track, 9729 cylinders
    Units = シリンダ数 of 16065 * 512 = 8225280 bytes
    
     デバイス ブート   始点      終点  ブロック   ID  システム
    /dev/hdc1             1       487   3911796   83  Linux
    /dev/hdc2           488       549    498015   82  Linux スワップ
    /dev/hdc3           550      5100  36555907+  83  Linux
    /dev/hdc4          5101      9729  37182442+  83  Linux
    

  • ファイルシステムの設定．

    $ sudo /sbin/mkfs.ext3 /dev/hdc1
    $ sudo /sbin/mkswap /dev/hdc2
    $ sudo /sbin/mkfs.ext3 /dev/hdc3
    $ sudo /sbin/mkfs.ext3 /dev/hdc4
    $ sudo /sbin/tune2fs -c 1000 /dev/hdc1
    $ sudo /sbin/tune2fs -c 1000 /dev/hdc3
    $ sudo /sbin/tune2fs -c 1000 /dev/hdc4
    

  • /etc/fstab

    /dev/hdc1               /mnt/hdc1               ext3    user,noauto     0 0
    /dev/hdc3               /mnt/hdc3               ext3    user,noauto     0 0
    /dev/hdc4               /mnt/hdc4               ext3    user,noauto     0 0
    

  • 次の命令は 絶対にやってはいけない もの

    $ sudo dd if=/dev/hda1 of=/dev/hdc1
    

    やってはいけない理由: partition table が壊されるから; やってしまったときの対策: 再 mkfs.ext3 しかない!
  • しょうがないので， sudo mount /mnt/hdc1 などしてから， まぬけっぽく

    cd /; for d in `ls -d /* | grep -v home | grep -v mnt | grep -v proc`; \
    do sudo mkdir /mnt/hdc1/$d; sudo rsync -av $d /mnt/hdc1/$d/ ; done
    

    そして cd /home && sudo rsync -av * /mnt/hdc3/ & …… あ， /home は dd わざ使ってもよかったんだ ……
  • うーむ， MBR コピーわざ とかあるみたいだな …… これはまた今度
• 作業途中に昼飯． 院生密度ほぼゼロ． さすがは大晦日．

[2006 年最後の日没]

本日の日の入りは 1609．
なんとも平穏のうちに．

• HDD にかまける作業はいったん中止して， 小山さんにサロベツ主要 8 植物種のやちパワー依存性を 推定する ZIP 的階層ベイズモデル関連のあれこれの納品 & モデルなどの解説． 今年の計算ぎょーむはこれにて終了 …… なのかな?
• ぢりぢりと主 HDD → 予備 HDD の差分バックアップ用の スクリプト書き． 以前にも書いたすごく便利 & 極めて危険な rsync まかせである …… 試験運転でさんざんしくじる． どーにかこーにかできた． /etc/cron.daily の「毎日の定期的な自動運転作業」のバックアップ用のやつを書き換える． 一件落着．
• このスクリプトにいちおう上記の MBR コピーわざ も含めておいた …… しかしながら， /dev/hda の MBR を /dev/hdc にコピーして動作するのかどうかは不明． まあ， 主 HDD がクラッシュしたときにその真価が問われるわけだが ……
• 1910 研究室発． あ， しまった． 今日は北 12 生協は 1800 店じまいか …… 地下鉄北 12 条駅セイコーマートまで戻ってみたり． 1940 帰宅． 晩飯．
• [今日の運動]
  • ついに年末は運動せぬまま ……
• [今日の食卓]
  • 朝 (0820): クロワッサン．
  • 昼 (1330): 研究室お茶部屋． 食パン．
  • 晩 (2050): 蕎麦． エビの天ぷら．