KuboLog 2007-09-05

2007 年 09 月 05 日 (水)

0730 起床．書類上の夏休み 3 日目，最終日．朝飯．コーヒー． 0915 自宅発．曇． 0930 研究室着．

伊庭さんから数セミ原稿についていろいろとコメントいただいたので，それにそった修正．おそるおそる編集部に連絡してみる．まあ，ダメだったら校正の段階であれこれ修正する，と．

どうも私の場合，「原稿みなおし」という作業に関してみなおし回数が増えるにつれて，劇的に粗雑化していくように思う．どうにかならんものかね ……

また ECOAS プランクトンデータの統計モデリングの検討．植物プランクトン体内での窒素 : リン = だいたい一定，を実現する方法がわからん …… たとえば，乱暴かもしれないけど，おおよそ No : Po = f : 1 ふきんをふらふらさせるにはにするためには

	dNo/dt = ... + r * Ni * No * (f * Po - No) + ...
	dPo/dt = ... + r * Pi * Po * (No - f * Po) + ...

とすればよいのかしらん?

それとも，そもそもこんなよけーなものつけなくてもいいのだろうか? …… いやいや，系全体が「平衡」状態でない場合であっても，植物プランクトン体内の窒素 : リンはめちゃくちゃにならないだろうから，何か上のような安定化のしくみが必要だよなあ．

とはいえ，すといきおめとりーだか何だか生物の体の中における「物質 A : 物質 B = 一定」なる考えかたってホントに正しいのか? あろめとりっくだとすると， log(A) : log(B) = ほぼ一定，ぐらいかもしれないじゃないか? ECOAS データの「観測値」の散布図でみるとそんなかんぢだ．だとすると，上の式は

	dNo/dt = ... + r * Ni * No * (f * log(Po) - log(No)) + ...
	dPo/dt = ... + r * Pi * Po * (log(No) - f * log(Po)) + ...

てなことになってしまうのかしらん? まあ，でもプランクトンというコンパートメントを考えると，やはり log ぬきのほうが妥当なのかなぁ．

しかし上のように窒素・リンの取りこみ速度をこういうふうに制約してしまうとぜんぜん動かなくなるような気がする．系が不安定ならいいんだけどね．やっぱり炭素動態も明示的に書かんといかんのかしらん? Chl-a の動態で代替できる?

	dNo/dt = ... + r * Ni * Ca * g1(Ca - fN * No) + ...
	dPo/dt = ... + r * Pi * Ca * g1(Ca - fP * Po) + ...
	dCa/dt = ... + g2(No, Po) * Ca + ...

この g1 だの g2 だのがややこしい関数になりそうだけど，まあどうしようもないよね．いやはや．

などと考えつつ昼飯．

とりあえず，水温データの階層ベイズモデル化．

> summary(as.matrix(post.mcmc)[, c("wt.a", "wt.d")]) wt.a wt.d Min. :-0.00477 Min. :-0.0127 1st Qu.:-0.00311 1st Qu.: 0.0443 Median :-0.00270 Median : 0.0591 Mean :-0.00270 Mean : 0.0586 3rd Qu.:-0.00222 3rd Qu.: 0.0744 Max. :-0.00119 Max. : 0.1147

[ただし，だ ……]

観測日があとになるほど標高の低いところに移動している，という交絡はあるんだよね．まあ，しょうがない．

あ，「年」変動いれるの忘れてた，と思って念のため調べてみると ……

> tapply(d$altitude, d$year, summary)
$`2003`
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
      5     115     250     357     660     845 
$`2004`
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
      5     232     660     687     852    2120 
$`2005`
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
    419     658    1130    1130    1530    2000 
$`2006`
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
    815    1450    1560    1670    1930    2700

…… 観測年によって調べている標高がぜんぜんちがう …… これが! 野外観測データだッ! (東北地方が中心になっているということで， BAOH ナレイションふうに叫んでください)

よくわからんまま，いちおう年変動もいれてみるか．あ，緯度いれるのも忘れてた …… ということでそのあたりモデルを変更してみた結果．

> summary(post.mat[, grep("^wt.[maldy]", colnames(post.mat))])
      wt.m           wt.a                wt.l             wt.d         
 Min.   :16.6   Min.   :-0.004162   Min.   :-0.704   Min.   :-0.00679  
 1st Qu.:17.3   1st Qu.:-0.002885   1st Qu.:-0.309   1st Qu.: 0.05031  
 Median :17.6   Median :-0.002494   Median :-0.197   Median : 0.07077  
 Mean   :17.6   Mean   :-0.002480   Mean   :-0.197   Mean   : 0.07051  
 3rd Qu.:17.8   3rd Qu.:-0.002044   3rd Qu.:-0.085   3rd Qu.: 0.09212  
 Max.   :18.9   Max.   :-0.000617   Max.   : 0.368   Max.   : 0.16420  
    wt.y[1]           wt.y[2]           wt.y[3]           wt.y[4]       
 Min.   :-2.4010   Min.   :-1.1920   Min.   :-2.0170   Min.   :-2.4560  
 1st Qu.:-0.2997   1st Qu.:-0.0318   1st Qu.:-0.0760   1st Qu.:-0.4274  
 Median :-0.0279   Median : 0.0800   Median : 0.0557   Median :-0.0720  
 Mean   :-0.1460   Mean   : 0.3046   Mean   : 0.1192   Mean   :-0.2371  
 3rd Qu.: 0.0890   3rd Qu.: 0.5906   3rd Qu.: 0.3653   3rd Qu.: 0.0377  
 Max.   : 1.7660   Max.   : 2.8530   Max.   : 2.1380   Max.   : 1.5140

[タテ軸は日・年補正した平均水温]

まあ，補正してもほとんどかわらん．変数名は wt.std より wt.adj とかにしたほうがいいかな?

…… というか，「補正」とかやらなくてもいいような気もしてきた．こういう山地の湖沼の生態系の回転速度がわからんのだけど，もし気温の日変化速度よりも窒素・リンの変動が速ければそれでいいわけで．うーむ，へいこう状態って何なのでしょう．実際には平衡でも何でもないんだろうなあ．

こういうやってもやらなくてもよさそうな水温のモデリングに手をつけたのは，他にも理由があって …… ゐんばぐす呼び出しのしくみを関数化によって多少はすっきりできないだろうか，というアイデアがあって，それをちょっと実装してみて試験してみたかった次第だ．下のごとく， library(R2WinBUGS) の wrapper 関数群を定義したら，多少はマシではなかろうか，と．

source("WinBUGS.R")
load("../data/d.RData")

set.data(list(
	N.site = nrow(d),
	N.year = 4,
	#
	W.temp = d$w.temp,
	Altitude = d$altitude - mean(d$altitude),
	Latitude = d$latitude - mean(d$latitude),
	Date = d$date - mean(d$date),
	Year = d$year - 2002, # 1:4
	#
	Tau.error.wt = 10,
	Tau.noninformative = 1.0e-3,
	Hyper.gamma = 1.0e-3
))

clear.param()
set.param("wt.m", 0)
set.param("wt.a", 0)
set.param("wt.l", 0)
set.param("wt.d", 0)
set.param("wt.y", rep(0, N.year))
set.param("d.wt", rnorm(N.site, 0, 0.01))
set.param("tau.wt.y", 1)
set.param("tau.d.wt", 1)
set.param("wt.adj", NA)

post.bugs <- call.bugs(n.iter = 5000, n.burnin = 1000, n.thin = 20)

これだってめんどくさいぢゃん，と思われるかもしれないけど，今までに比べれば「圧倒的にすっきり」 (当社比) 書けてるのである．しかも，この見かけより重要なのはモデリング試行錯誤がやりやすく (当社比) なってる，というあたりだ．

1815 研究室発． 1830 帰宅．体重 68.0 kg．プランクトンやせ，か．うんどう．晩飯の準備．晩飯．

ベイズモデルでどうやればすっきりモデル選択みたいなことができるのか，いまいちよくわからない．そういう研究はいろいろとある．たとえば， mixture prior を使う，といった方策だ …… これって実装が少しばかりめんどくさそうだ．

[今日の運動]

エアロバイク 45 分間． 400kcal．

[今日の食卓]

朝 (0800): 米麦 0.5 合．ニラ・ニンジン・ショウガ・エノキダケ・豚肉のカレー．
昼 (1340): 研究室お茶部屋．ヨーグルト．リンゴ．
晩 (2100): 米麦 0.8 合．ネギ・ピーマン・豆腐・卵の炒めもの．キャベツ・ネギのサラダ．コマツナのゴマあえ．キュウリ．

ぎょーむ日誌 2007-09-05

2007 年 09 月 05 日 (水)