ぎょーむ日誌 2006-11-(01-10)
2006 年 11 月 01 日 (水)
-
0810 起床．
なんかふちょー．
0830 自宅発．
晴．
0845 研究室着．
朝飯．
コーヒー．
-
また地衣類計算やらせつつ，
(昨日かいたことのふぉろーで)
R のインストール
というのをごくごく簡単にまとめてみたり．
-
来週火曜日が提出しめきりの査読ぎょーむがあるので着手するか．
今回はいつもよりはまともそうなやつなんだけど，
そのぶん読むのがたいへんかもしれない．
ちょっと物理学よりなやつだし．
-
すぐに脱線してメイルかきとか．
-
昨日からぼちぼちと続けている
WinBUGS の
n.iter = 60000, n.burnin = 30000
トリプルチェイン計算のまとめ
(一回の計算時間は 3000-4000 秒ぐらい):
| random effects |
「明るさ」 |
DIC |
めも |
| 個体 / 岩面 / 岩 |
(NA) |
4007.0 |
|
| 個体 / 岩面 / 岩 |
logit |
4006.1 |
|
| 個体 / 岩面 |
(NA) |
4003.3 |
|
| 個体 / 岩面 |
logit |
4001.5 |
このへんか? |
| 個体 / 岩面 |
log |
4001.8 |
収束わるい |
なぜ「岩ぬき」が良いかというと，
この志水さんの観測データではたいていの場合
「ひとつの岩から岩面 2 枚」
というサンプリングになっていて，
しかもその二枚の岩面はかなり異なるので
「岩ごとの効果」
の事後分布が平均ゼロでありながらばらつき大きい，
という状況にあるため．
別のいいかたをすれば，
岩面ごとの「個性」がきついのに，
岩面数が 2 では「岩ごとの特徴」を推定しづらい，
ということ．
-
R2WinBUGS
でもって Linux 上のゐんどーづアプリケイションソフトウェア実行環境
Wine
上で WinBUGS よびだして計算させると，
MCMC 計算終了時にこういう図をだして止まる
(
bugs(..., debug = TRUE) 指定してる場合)．
階層ベイズモデルの MCMC 計算の 3 反復チェインの遷移の一部，
というところか．
これは上の表で「収束わるい」と書かれている log(明るさ)
なモデルの計算結果である．
収束よさそうに見える上のふたつは岩面 62 番と 63 番の「個性」
(random effects 事後分布)
……
63 枚ある岩面の中で 62 番は
(random effects なので理由はまったく不明ながら)
地衣類の繁殖には特に不適そうな岩面，
ということをあらわしている．
-
で，
bugs()
から受けとった bugs オブジェクトを使って
R の中で図を描いたり，
あるいはその bugs オブジェクトを save()
したり，
と．
-
もはや WinBUGS の計算もだいたい上のように終わりつつあるので
……
そうだな，
同じ条件で JAGS で計算やってみるか．
-
JAGS だと 53 分かかった
……
ちょっとサンプル数指定をまちがったんだけど，
JAGS で 90K MCMC step の計算をやったので
10K MCMC step あたりだと 350 秒ぐらい．
いっぽう R2WinBUGS の WinBUGS だと 180K step の 3500 秒ぐらいだから，
10K MCMC step あたりだと 190 秒ぐらい．
なぜか WinBUGS のほうが速い
(アヤしげな魔術?)．
混交のぐあいに関しては今回はどちらも同じようなかんじ，
か
……
いや，
なぜかしら WinBUGS のほうがしゃかしゃか混ざってるようなかんぢだな
(アヤしげな魔術?)．
-
昼飯．
-
地衣類に関する MCMC 計算はひとまず終わったような気がするので，
rsync
してから，
計算につかってた Dell 機をネット経由で shutdown -h now．
ここから先はまた志水さんとの相談が必要になるわけだが
……
ここ二週間ぐらいお会いしてないので
(他の仕事がおいそがしいのでしょう)，
これ以上はススめようがない，
と．
-
あ，
例のドイツ語版からの翻訳本
「Rの基礎とプログラミング技法」
ってすでに
入手可能
なのか
……
いやはや，
すでに「在庫切れ」ぢゃん．
-
査読ぎょーむつづき．
ざっとながめてみたら，
この方式の計算の意図の概要はつかめた
……
というのも，
これって私もさんざん考えていたんで．
で，
さんざん考えたあげくに
当方が使わなかった理由は観測データと対応がつかないから．
今回みている原稿はそのあたりは斟酌していない
……
まあ，
それは「趣味のちがい」と気にしないことにして，
だ．
なンかこの方式の計算方法にはすごい
「読みぬけ」
があるような気がするんだよね．
それが何だかわからんのだけど．
-
一点ヘンなところというか，
説明が苦しくなるところに気づいた．
(こういう問題考えたことあるヒトならすぐに気づく)
めんどうを回避するべくうまく条件つき確率を定義しているわけだが，
それゆえに
「えー? なんでこれは……」
なる素朴な疑問に答えられない，
と．
「だって他のヒトたちもそうしてるもん」
はありうるが
……
いや「他のヒトたち」の計算方式とは少し別ものだと
今さらながら気づいた．
これは興味ぶかいな．
指摘する価値はある．
-
まあ，
作文なんかは私なんかよりもうまいわけで，
このへん指摘したところで
``good question''
とか何とかでてきとうにごまかされそうな気もしますが．
作文能力が低い私だったら，
そういう指摘されたら計算のやりかたのほうを変えるけど
(そっちのほうがラクなので)．
-
いずれにせよ，
かなり抽象化された「あっちの世界」のハナシなので，
どういう現実ばなれも許されるのである．
つまりそれにあれこれ指摘してみるってのは
……
いわば，
ドラえもんの道具が物理法則から逸脱してるように見えるという理由で，
作者に「修正しろ」と要請してるようなもの，
かしらん．
うーむ，
すると
「そこで『どこでもドア』を使うのはおもしろくないので，
むしろ『タケコプター』にするといいハナシになりますよ」
ぐらいしか言えないではないか?
-
などと現実逃避のらくがきしてるヒマがあったら，
この長ったらしい原稿のすみずみまで熟読すべきだよな．
たしかに作文はうまいし，
巧妙に難点を回避しているけれど
(というか回避してなかったらそもそも計算できんわけだが)，
その回避動作そのものに上のごとく何やら墓穴系なところがある．
-
……
と書いてみてから，
もういっぺんだけこのヘンな点を検討して
……
みると，
「この世界のルール」内であっても
(「どこでもドア」より
「タケコプター」のほうが移動速度が速いかのような描写がある，
といったたぐいの)
「矛盾」をきたしている可能性があることに気づいた
……
そして，
なぜ気づいたのかといえば，
このあたりもまた大学院生のころにまったく同じことを考えていたからだ，
とこれまた今になって思い出した．
-
まあ，
この著者はそのあたりなかなかヌケ目がないヒトみたいなんで
(というか計算やった本人ならとうぜん気づいてるよな)，
この長い作文のどこかで
回避策・正当化・いいわけ
(「じつは『どこでもドア』
は任意の二地点を自由に接続できるわけではなく……」
といったたぐいの)
がちゃんと書いてあるかもしれん
(このへんの回避策らしきモノは私も考えたことがある)
……
ということで，
おとなしく虚心坦懐にぢりぢりと読みますかね．
しかしその前に「考える」ことで
(というより「思い出す」というべきか?)，
このあたりの「土地勘」がよみがえってきたかんぢではある．
-
夕方から，
お茶部屋にて松田さんと駒ヶ岳気象データ変換作業．
ちょうど一年前
にも同じことやっているな
……
ぎょーむ日誌の簡単な記録であっても，
何もないよりはマシ．
-
そして昨年つくった Perl フィルターが意外とよくできていて，
今年の欠側値だらけのデータロガー出力にも対応できているかのように見える
……
最後は R
の
read.csv()
で読ませるわけだが，
この関数もなかなかよくできているというべきか，
「列が少ない行」は自動的に NA
でうめているな．
-
1800 あたりをすぎると，
お茶部屋は大学院生
(これは，
より具体的には，
太古よりこのお茶部屋を支配してきた特定性別がわのヒトたちなんだけど，
このへんをうかつに書くとキビしく叱責されるので……)
の密度が高まると異様な喧騒につつまれるかんぢだ．
「食べものに興奮する大学院生」
といったことをこのへんに書いたりしたんだけど，
よく考えてみるとこれは単純化しすぎで，
おそらく
「大学院生は他の大学院生に興奮している」
(密度依存性が重要)
という部分も含めたほうがより正確な現象の記述になっているのかも．
-
フィルターを何ヵ所か修正して，
駒ヶ岳気象データひとまずは終了．
1940 研究室発．
2000 帰宅．
晩飯．
-
2321 地震．
震度 3 ぐらいか．
いやそんなに
ゆれてないか．
-
[今日の運動]
-
[今日の食卓]
- 朝 (0850):
研究室お茶部屋．
セイコーマートのサンドイッチ．
- 昼 (1330):
研究室お茶部屋．
食パン．
- 晩 (2200):
米麦 0.8 合．
ダイコン・ニンジン・ネギ・豆腐・煮干の味噌汁．
2006 年 11 月 02 日 (木)
-
0720 起床．
朝飯．
コーヒー．
0900 自宅発．
晴．
0915 研究室着．
-
駒ヶ岳気象データ変換，
やはり再履修通知が
……
まあ，
よくあることだ．
で，
もとデータをよくみると，
昨年度変換したデータもまた変わっていて，
観測時期によって
「コンマくぎり」「タブくぎり」
データが混在している!
とわかった．
おそらくこれは私を鍛えるために大学院生が設定してくれた試練なのだろう．
-
そして，
この観測時期によって異なる欠側データ，
というか観測時期によって測定してたり測定してなかったり，
という面倒をやっつける Perl フィルターが必要とされてるわけで
……
-
と苦闘してると，
とある事情の昼飯招集がかかり，
ぞろぞろと (百周年記念会館 1F) きゃら亭にいって昼飯．
-
昼飯後もフィルター工夫．
あまりトリッキーなわざで面倒を片づけるのは好きではないんだけど，
ややトリッキーに解決．
一年後にみたら理解できんだろうな
……
ということで自分のために簡単に注釈しておくと，
長期欠測しうるデータ列は全部「右側」によせてしまう，
という方法だ．
-
とりあえず標高ごと全データの図を描かせる
R
プログラムをかいて，
ざっとチェックしたことにする
(下の図はとある標高のもの)
……
まあ，
いいでしょ．
ということで，
松田さんにひきわたし作業．
-
査読ぎょーむのつづき．
昨日よみこめていなかったところはだいたい予想どーりだった．
しかしながら，
今回の査読報告かきは時間かかるかも
……
もう原稿みなくても，
すべての計算内容を自分の手で全部再現できるほどアタマに入ってる，
にもかかわらずだ．
-
計算方法はそうとうにヘンであるにもかかわらず，
この計算そのものは「正しい」
……
ただし，
この「正しい」の意味は
「文章で説明してることと数式にだいたい矛盾ない，
その説明は計算機プログラムとしてゴマカシなく実装可能，
計算結果あれこれの統計学的パターンも再現可能だろう」
ぐらいの意味で
(私のところにまわってくるような数理モデル系ダメ論文なんかは
このあたりもアヤしいのが多い)，
その結果の解釈にも資意的なところは少ない
(これまたヒドいのがまわってくることが多いんで)．
-
主著者は物理学系のヒトでそちら方面の趣味が露出したものだが
……
まあ，
この点はそもそもヘンなぢゃーなるからの査読依頼なので，
問題ナシ．
共著者の生物学者が書いたらしき引用あれこれは，
ややマトはずれなものだが
……
まあ，
これはよけーな作文を削れば改善されるだけのことだ．
-
気になる点がいくつかあって，
相互に連関している．
これをうまく査読報告として作文するのは面倒だからイヤだなあ，
とうとましく思っているわけだが
……
気になる点のひとつは著者が計算方法の
「ヘンさ」
をどこまで理解してるのか
よくわからないこと．
このヒトは勉強ずきらしく，
何か計算プロセスを説明するたびに
「これは以前の誰それの計算と同じ」
と書きたがるわけだが，
じつは核心部は解釈不可能なほどに，
何というか独創的にヘンなのである．
それをわかってないのか，
それとも先方も面倒だから「前と同じでしょ」と書いてるだけなのか．
このへんは昨日かいてた回避→墓穴問題まわりと直結している．
-
そして
「従来方式の発展で今回の計算を」
とも言いたがるわけだが，
この従来方式の
「踏襲する必要がないところまで踏襲した」
ために best な方法になっていない部分がある，
とわかった．
これは「カタチだけマネした」つまり考えが足りないせいなのか，
新しい定式化のアイデアがないだけだったり，
「そりゃわかってるけど結果さえだいたい都合良ければよかろうなのだ」
と開きなおってるいるのか
……
えーい，
もうちょい計算方法がマシになるよう努力してくれていたら，
こういう面倒な指摘はやらなくてすんだのに．
-
1830 ごろからお茶部屋の小パーティーに参加．
主賓は IGBP
GLP
札幌ぶらんちの新事務局長さん
(ナイジェリア出身の人)，
かしらん．
-
終了後も上の査読問題がアタマから離れず，
夜の A 棟 8F 廊下をうろうろと巡回．
-
やはりあの計算まわりに関する理解はまちがってないよなあ，
と思いつつも本日は
「計算は正しい，
しかしヘンだということを十分説明してないし，
best でない部分がある」
といっためんどう作文する気にもなれず撤退．
2055 研究室発．
2110 帰宅．
-
2310 自宅発北大構内走．
晴れてるせいか寒い．
そしてアタマの中は査読計算にヤラれたまま．
2355 帰宅．
-
[今日の運動]
-
[今日の食卓]
- 朝 (0830):
米麦 0.6 合．
ダイコン・ニンジン・ネギ・豆腐・煮干の味噌汁．
- 昼 (1220):
きゃら亭．
A ランチのサケのムニエル．
700 円．
サラダ・コーヒーつき．
- 晩 (1840):
お茶部屋なべもの．
2006 年 11 月 03 日 (金)
-
0850 起床．
朝飯．
コーヒー．
怠業．
洗濯．
怠業しつつ査読報告の構成を整理してみる．
1050 自宅発．
曇．
1105 研究室着．
-
査読ぎょーむつづき．
査読報告かきに着手．
-
しかしすぐに脱落して，
Vine Linux 3.2 の
sudo apt-get update && sudo apt-get dist-upgrade
などやってみると
……
お，
谷村さん packaging の R-2.4.0
がインストールされた．
-
で R でも
update.packages(ask = F)
すると
……
library(VGAM)
もびみょーに更新されており，
これでついつい遊んでしまう．
[lvplot()]
library(VGAM)．
オランダのクモデータ
(
data(hspider))
の constrained additive ordination (CAO) model
を
cao() で推定させたのを図にするとこうなる．
-
で，
査読報告を完成させようとしたんだが
……
-
完成にはほど遠い状態で敗退．
1925 研究室発．
1940 帰宅．
体重 72.6kg．
晩飯．
-
で，
明日は
徳舜瞥
(とくしゅんべつ)
登山のおとも．
-
[今日の運動]
-
[今日の食卓]
- 朝 (0900):
米麦 0.7 合．
ダイコン・ニンジン・ネギ・豆腐・煮干の味噌汁．
- 昼 (1520):
研究室お茶部屋．
食パン．
- 晩 (2010):
うどん．
ダイコン・ニンジン・ネギ・豆腐・煮干の味噌汁．
2006 年 11 月 04 日 (土)
-
0550 起床．
0600 自宅発．
曇．
0615 研究室着．
朝飯．
うすいコーヒー．
-
0705 また甲山さん・川合さん山行のおともとして地環研前出発．
本日は
徳舜瞥
(とくしゅんべつ)
山．
かーなびのご指示どーりの，
中山峠ごえの喜茂別で支笏湖方面へ反転の大滝村
(行ってみたら伊達市に合併されてた)
いきコース．
[08:39 羊蹄山]
見えにくいけど画像中央あたり．
喜茂別町にはいったあたり．
この奇妙な山名の由来は
Wikipedia
など参照．
とにかくこのへんでは突出して高い山 (1898m)．
このあと，
(細い林道はかーなびの世界では「なかったこと」にする)
かーなびにだまされて登るべき林道をまちがえる．
[09:48 登山口からみる徳舜瞥山]
正しくは紫明川ぞいの林道を標高 730m までのぼる．
0930 登山口着．
徳舜瞥は大滝周辺では突出して高い山．
まあ，
この「うしろ」 (南方) に位置するホロホロ山のほうがびみょーに高いけど．
本日は紫明川ぞいの林道から徳舜瞥北尾根にとりつき登るルート．
[09:55 北尾根の谷側を登る]
札幌あたりからの登山客が多いためか，
ルートはよく整備されている．
[10:27 七合目, 標高約 950m]
稜線上にあがる．
ちなみに登山口は五合目 (標高約 730m)
とのこと．
[11:10 九合目, 標高約 1150m]
ハイマツ帯から旧大滝村を見おろす．
天気は晴れてて，
羊蹄山・尻別岳・支笏湖・恵庭岳などがみえる．
この北尾根の一番下には大滝の「飛行場」がある
……
といっても軽飛行機専用の小さいやつだけど．
[11:21 徳舜瞥山頂付近]
むこう (東側) に見えてるのはホロホロ山．
尾根上のルートは特に急登もなくなんとなく歩いてるうちに登れてしまった．
(
やや大きい画像)
[11:22 徳舜瞥山頂，標高 1309m]
雨ふりの空沼岳，
雪にみまわれた札幌岳とは異なり，
今日は好天にめぐまれた．
気温もなぜか高い．
雲はそれなりにあったので，
羊蹄山・恵庭岳などは見えにくかったけれど．
1130 徳舜瞥山頂発，
すぐ東にあるホロホロ山めざす．
[12:00 ホロホロ山への縦走]
ホロホロの山頂からみた縦走路．
うしろ (西側) に見えてるのは徳舜瞥山．
この縦走路は北海道には珍しくちょっと岩っぽい尾根．
(
やや大きい画像)
[12:01 ホロホロ山頂，標高 1322m]
ホロホロ山は山頂直下が急登になっていて，
「登った」という実感がある．
登山者のみなさん，
たいていはこちらまで縦走するようで人口密度たかかった．
ここで昼飯くってから，
もときたコースをひきかえして下山．
[14:02 駐車場に下山]
……
ということで 4 時間と 20 分ほどで登山終了．
冬山状態になるんでは，
という懸念とはうらはらに「お手軽な秋山」とういかんぢでした．
1410 駐車場発．
[14:30 本日の温泉]
蟠渓 (ばんけい) ふれあいセンターなる
壮瞥町の施設．
あまり大きくない温泉．
リンゴ (レッドゴールド) 12-3 個 500 円だったので思わず購入．
1515 同発．
……
によった．
きのこ・野菜の直販 & 軽食など．
大学院生の食欲に圧倒される．
甲山さんはキノコつめあわせなど買っておられる．
国道 276 号ぞいにある．
隣接する「きのこ村」と競合しており，
さらに
裁判ざた
にもなってるとか．
[16:41 支笏湖と恵庭岳]
帰りは中山峠はとおらず
(定山渓あたりが混雑するのでは，
という前回の札幌岳帰路の戦訓)，
支笏湖まわりのルートから帰る．
恵庭岳もよくめだつ突出して高い山なんだが
……
標高は 1320m で本日のぼった徳舜瞥山・ホロホロ山とそんなに変わらない．
-
札幌市内の創成川ふきんはやや渋滞．
また甲山さん車で大学まで送っていただいて，
1755 北大着．
本日もありがとうございました．
-
そしてお茶部屋では今年最後 (?) の苫小牧調査その他を片づけた
M1 の皆さんがハイになっている．
-
うだうだと山行ぎょーむ日誌なんぞを書いてから，
2050 研究室発．
2105 帰宅．
晩飯．
-
[今日の運動]
-
[今日の食卓]
- 朝 (0625):
研究室お茶部屋．
食パン．
- 昼 (1205):
ホロホロ山山頂．
せぶんいれぶんのにぎり飯．
今回も梅．
- 晩 (2130):
スパゲッティー．
レトルトパウチドのナス・トマトソース．
2006 年 11 月 05 日 (日)
-
0720 起床．
雨．
朝飯．
コーヒー．
怠業．
洗濯．
-
1255 自宅発．
曇．
1310 研究室着．
-
夜中までデータ解析こんさる．
-
まずは，
R
によるアヤしげ数値積分のやりかたの質問．
こういう計算できる関数つくって返信．
-
1400 から夕方まで農学部の高橋さんととある発芽データ解析．
単純な logistic 回帰だが
……
plot 間の「個体差」をくみこんだ
GLMM
にすると統計モデルが格段によくなったので感銘うけた．
-
作者 Göran Broström さんと熱心なユーザーたちによって
library(glmmML)
はぢりぢりと改良されており，
-
family = binomial において
formula = cbind(生起した回数, 生起しなかった回数)
~ ...
書法が使えるようになった
(ふつうの glm() と同じように書ける)
-
stepAIC()
を使った自動モデル選択ができる
-
「切片」の最大事後確率 (MAP) 推定値も計算される
……
といったことが最近はできるようになっている．
データ解析ひとまず終了後，
glmmML ペイジ
をちょっと作ってみる．
-
しかし曜日に関係なく，
お茶部屋では院生が常に食いもののハナシばかりやってるわけで
……
-
夕方から夜にかけて松田さんのデータおてつだい．
例によって例のごとくと申しますか，
ゑくせるゆーざーにありがちな「横長」経時データ構造 128 列
……
こいつをまっとうで使いやすいタテ型に変換せよ，
という問題．
-
1800 個体ぶんのヨコ→タテ変換じたいはまったく簡単なのだが
……
簡単でないのは R
における処理時間の短縮である．
なまなかな方法では 20 分かけても終了しない．
そこで私がいろいろと工夫して
sapply()
三段がさねの超長 vector 生成おりかえして matrix 化，
というふつーならざるアラわざを使っても 350 秒を費した．
-
ヨコをタテにするだけなのに R ではなぜこんなに時間がかかるのか?
理由はおそらく
「R 内部における
data.frame()
の動的なメモリ確保」
といった問題に関連していると思われる．
データ量が増えると「加速度的に」処理時間が増大し，
10 個体ぶんなら 1 秒で終わるから，
100 個体なら 10 秒で終わるかとおもいきや 12 秒かかり，
ならば 200 個体ならとためしてみると 35 秒かかったりする．
変換もともしくは変換さきの
data.frame()
が肥大化してくると一処理ごとに費される時間がすごく増大する，
という構造がある
……
-
このあたり解決できぬまま撤退．
2245 研究室発．
雨．
2300 帰宅．
晩飯．
-
やっぱ，
ケガれ言語 Perl はこういう問題に強いわ
……
(あまり公平な比較ではないけれど)
R でやると 350 秒を費す変換，
Perl でやると 2.4 秒以下で終了．
とりあえず，
変換問題じたいは解決．
-
ここで R vs Perl の変換スクリプト比較をしてみると
……
おそらく R においても
「変換したらすぐに吐きだす」
という方式にすれば 10 倍ぐらいの高速化は実現するのではないか?
-
そして夜中のプログラミングは健康にわるい．
-
[今日の運動]
-
[今日の食卓]
- 朝 (0800):
バターロール．
- 昼 (1215):
明太子スパゲッティー．
- 晩 (2350):
米麦 0.7 合．
ネギ・ピーマン・卵の炒めもの．
2006 年 11 月 06 日 (月)
-
0800 起床．
朝飯．
コーヒー．
0850 自宅発．
雨．
0905 研究室着．
-
昨晩のデータ変換問題にまたしつこくとりくむ．
R
だと 350 秒かかっていた変換処理時間を
5 秒
程度に短縮できた
(Perl の 2.4 秒にはまだ及ばないが)
……
コツは
「いっさい
data.frame() を使わない」
こと．
すなわち readLines()
で読んで変換してただちに
writeLines()
ですぐに吐きだす，
という方式．
-
文字列変換というオマケをつけると 14-16 秒ぐらいまで増大してしまった．
しかも R ではうまく扱えない文字コードが
……
try(gsub(...))
でごまかす．
-
writeLines()
はいろいろ不便なので
out <- file(file.output, "w")
...
cat(..., file = out)
...
close(out)
と cat()
で吐きだす方法にした．
-
しかし松田さんのところで作業してるうちに
gsub(...)
による置換が増えてしまった，
これがえらく時間がかかるので
(R の弱点は文字列処理か?
……
後記:
あとになってようやく気づいたんだが，
data.frame 一行ずつ gsub(...)
やらせるのは最悪な処理だ，
というのも行数ぶん正規表現を「こんぱいる」することになるから;
教訓:
gsub(...) などは「全データかためて」ほどこすべし)
……
けっきょく松田さん PC に
ActivePerl
をインストールして Perl で処理することに．
多少のめんどうな置換 (コンマ濫用対策とか)
をやらせてもなお，
やはり Perl は「むちゃくちゃ」速い．
-
おそい昼飯．
井田君と R
の
data.frame()
のハナシ．
たとえば列を置換する場合には，
sapply()
で「いっぱつ置換」しないと使わないと遅くて
(あまり良い例ではないかもしれないけど)
df$x <- sapply(df$x, function(v) ifelse(v < 50, "A", "B"))
といったふうに書かないといけない．
-
ついでに FAQ への回答．
「R で factor の水準の順番を変える
にはどうしたらよいか?」，
これは簡単で以下のごとし．
> (treatment <- factor(c("1", "2", "control")))
[1] 1 2 control
Levels: 1 2 control
> factor(treatment, levels = c("control", "1", "2"))
[1] 1 2 control
Levels: control 1 2
あるいは最初に factor(...)
するところで水準を指定すれば簡単だが，
> (treatment <- factor(c("1", "2", "control"), levels = c("control", "1", "2")))
[1] 1 2 control
Levels: control 1 2
このワザを使いたくなるのは read.csv()
で読んだデータの水準変換だからね．
その場合は上とおなじくいっぱつ置換．
> df$x <- factor(df$x, levels = c("control", "1", "2"))
-
そして昼めし後も松田さんデータ変換を少々
……
-
うーむ，
「あの」
Ripley
先生が東京に来る!
12/8 (金) & 12/9 (土)
の統計数理研究所共同研究
「Rの整備と利用」
研究会
……
まあ，
土曜日のほうは私も
MCMC
計算の苦闘についてなんぞをしゃべる仕儀におちいったわけだが
(ということで，
先月後半は R2WinBUGS
などを研究していた次第)．
-
さてさて，
12/8 (金) の内容かなり計算機統計学よりで難しそうだが
……
やはり Ripley 御大の講演も聞いてみたい気が．
となると 12/7 (木) に札幌発の 12/10 (日) に東京発といった日程になるな
……
だいじょーぶなのか?
-
うう，
昨日からアタマが「下うけモード」になってしまって，
査読報告かきに復帰できん
(と他人のせいにしてみる)
……
-
けっきょく査読報告かきは進捗せぬまま撤退．
お茶部屋雑談にとらっぷされる
……
するとマレイシアで調査中の矢澤さんからお電話が．
要点のメモをとり，
のちほど甲山さんにメイルで報告，
と．
2025 研究室発．
雨．
2040 帰宅．
晩飯．
-
[今日の運動]
-
[今日の食卓]
- 朝 (0820):
米麦 0.7 合．
ネギ・ピーマン・卵の炒めもの．
ワカメスープ．
- 昼 (1340):
研究室お茶部屋．
Sanks にぎり飯．
- 晩 (2200):
米麦 0.7 合．
チンゲンサイ・ハクサイ・豆腐のカレー．
2006 年 11 月 07 日 (火)
2006 年 11 月 08 日 (水)
-
0700 起床．
ねむい．
朝飯．
コーヒー．
0850 自宅発．
曇．
0905 研究室着．
-
昨日まで放置してたメイルかき & ネット雑用，
など．
-
マレイシア続報．
なぜたくさんの小径木が死んでいたかという理由が判明．
甲山さんが「1999 年のデータ」と言ってたファイルはじつは
1995 年に観測されたものだった!
そりゃ 11 年 (以上) もたてば小個体ならそれなりに死ぬでしょうよ
……
教訓: Pasoh にかぎらず他人から渡された毎木データは
「いつどのように取られたものか」
をしつこくしつこく確認すること．
-
粕谷さんから多変量解析のアヤしさについて考察せよ，
という宿題がでてたので，
そのあたりを考えてみる
……
というか PCA だの CA だのそもそも計算方法がアタマから
脱落していたのでそのあたりから復習．
この計算じたいはどちらも単なる Lagrange の
未定乗数法
で何も難しくはない．
-
私が提出したごくごくいいかげんな多変量雑談．
はい，固有ヴェクトルの成分の確率分布などはあまり考えられていませんね．
同じ研究室の露崎さん (後述する植生学で *CA のたぐいをよく使っている)
から固有値の寄予率などは「検定」する (おそらく bootstrap 的な手法でしょ
う) ことが近ごろでは求められている，といったハナシをきいたことがありま
す．固有ヴェクトルの成分のばらつきまでは調べてないような気がしますね．
...
「もともと」はピアソン (親父) あたりが始めた多変量正規分布の分散・共分
散行列の推定が起源でしょうから(安藤洋美「多変量解析の歴史」とか)，母集
団みたいなものは考えていたと思います．もしデータが多変量正規分布からの
ランダムサンプルであり，そこから主成分分析 (PCA) によって得られた固有
値・固有ヴェクトルのあつまり (行列) をひねくれば分散・共分散行列になお
すことができるだろうと思うのですが，そのやりかたはよくわかりません．
Lagrange の未定乗数法をもちいる現在の *CA のたぐいの起源は 1950 年代ぐ
らいではないかと思うのですが……まあ，こういう手法が開発された理由は，
分散・共分散行列を推定してながめてもよくわからないので，固有値・固有ヴェ
クトルを使って「主成分」を示すほうがわかりやすい，と考えられたためなの
でしょう．
勝手な憶測をつづけますと，計算方法から多変量正規分布という概念がぬけて
しまった時点で，粕谷さん指摘されるように「母集団を想定しない」不思議な
統計学的手法になったのでしょうかね．
あるいは三中さんが以前に言われてたように，このへんは情報圧縮 (でしたか?)
みたいな技術として使われてる，と．
最後に「例」の報告ですが，粕谷さんが探しておられる「ヘンな例」の範囲か
らも逸脱してるような気がするのですが，「植生学」は *CA な手法にかなり
依存している分野としてあげることができます．これはいろいろな場所で，わ
くをとってそのわくの中の植物種ごとの「植被度」なる「割合みたいな量」を
要素とするヴェクトルをつくり，それを *CA します．
おもしろいのは，きわめて植生学的には理想的な (本来なら「第一主成分」だ
けで完璧に説明できるはずの) データであっても，PCA，CCA などはありもし
ない「第二主成分」が推定される，ということです．図としては
http://hosho.ees.hokudai.ac.jp/~kubo/log/2004/1201.html#vegan
あたりに示しています．これは「arch 現象」とか呼ばれてて植生学のヒトた
ちのあいだでもよく知られています (そこで DCA なる想像を絶した奇怪な手
法が「発明」されてしまうのですが，それはまた別のハナシ)．
R
の VGAM
package などはこのあたりの複雑怪奇さにイヤけがさし，
それに対抗するべく
「多変量解析の統計モデル化」
をねらったものだと思うんだけど
(先日の lvplot()
作図例も library(VGAM))
……
じつは私はこのあたりはまだまだ不勉強な状態だ．
さてさて，
これから生態学まわりでの多変量解析はどうなっていくのであろうか．
てなこと書いてたら，
私が植生学的多変量解析をアヤしげよばわりしていると
(おそらく霊感か何かで)
察知した院生が，
私は自分の植生データを DCA で解析するつもりです，
といきなり宣言していった
……
まあ，
今日はダメでも
「いつかは」
このあたりきちんとした統計モデル化の勉強をやろう．
要するに「相互に独立ではないかもしれぬ成分たち」を要素とする vector の
「ぐるーぷ化」だの時間変化のモデリングだろう?
階層ベイズモデルでけりがつく問題だ．
私が目次 mail 配信うけとってる数少ないぢゃーなるのひとつ
Environmental and Ecological Statistics
(目次を mail で受けとらないとその存在すら忘れてしまいそうなので)，
今回はちょっと役にたちそうな情報もあり．
A spatial zero-inflated poisson regression model for oak regeneration
Author(s): Stephen L. Rathbun, Songlin Fei
Page: 409 - 426
DOI URL: http://dx.doi.org/10.1007/s10651-006-0020-x
修論とかの下うけ計算で使うことになりそう，
と思ってダウンロードしてみたら
……
やはりというか，
階層ベイズモデル + Gaussian Random Field
(これを spatial probit model として
「適・不適」隠れ値生成に使っている!)
なモデルの MCMC 計算 + EM アルゴリズムあわせわざ推定，
でした．
ということで，
皆さんあきらめてこのあたりも勉強されてください．
昼飯．
昼飯後はなぜか社会性昆虫のあれこれなどダウンロードして勉強．
1500 よりここ A 棟 5F 509 室で
社会性昆虫勉強会
に参加．
[社会性昆虫勉強会（社昆勉）札幌大会]
日時：2006年11月８日（水）15時より（懇親会１９時〜）
会場：北海道大学・大学院地球環境科学研究院 A509室
演題：
前川清人（富山大・理）
「シロアリのカースト分化と器官形成に関する研究」
柴尾晴信（東大・総合文化）
「社会性アブラムシにおける階級分化と社会制御」
-
主催者の
三浦さんの楽しげなごあいさつと，
その師匠である松本さんによる
(いままでは東京で開かれてた)
社昆勉 30 年間の略歴が紹介されたあと
(桑原万寿太郎，
なんてひさしぶりに聞いた)，
本日の講演はじまる．
-
で，
3 時間ちかくにわたって
前川さん
と
柴尾さん
のたいへん内容もりだくさんのハナシがあったわけだが
……
[びっしりメモ]
いやいや，
ちゃんと内容はフォローしてました，
と
……
前川さんの研究はシロアリのワーカー・ソルジャーへの分化によって
何が生じるのか (個体内)，
どのようにそうなるか (個体間?)．
柴尾さんは真社会性アブラムシの分化誘導抑制 (個体間相互作用)
- 齢差分業 (「お年より」に危険な仕事を!) - 化学信号への齢差応答，
といったハナシの流れ．
-
お二人のハナシや一連の質疑応答からはちょっとハズれるかもしれぬが，
拝聴してるうちにいくつか励起されてしまった，
例の東さん科研がらみのしみゅれいしょん研究だか何だかと関連しそうな
妄想みたいなものを整理してみよう，
と
……
-
sibling 相互作用．
シロアリでもアブラムシでも
「兵隊は (幼虫・ワーカーの兵隊誘導による) 兵隊増産を抑制する」
sibling 相互作用は普遍的というか，
「まあそうだろな」的に見られる．
では兵隊誘導はどういった個体間相互作用によるものか?
これも sibling 間相互作用だとするとすっきりするので，
そう考えたいわけだが
……
-
親による誘導，
というハナシもあるようで．
親が「こいつは兵隊，こいつはワーカー」
とか識別できてるのかしらん?
-
さらに，
教科書的な表現的多型ばかりがすべてを決めてるわけでもなく，
そもそも「何に分化するか」は遺伝学的背景も関与するらしい，
といった指摘があったり．
-
シロアリコロニー創設実験をやると，
創設期には兵隊は一匹だけ出てとりあえずそれ以上は増産されない．
あとで前川さんにお尋ねすると，
とりあえず一匹いればコロニー防衛
(シロアリの場合は他コロニーのシロアリ排除が重要だそうで)
はできるのでしょう，
と．
そるじゃーってワーカーよりも格段に「戦闘能力」が高いのか?
-
そるじゃーとわーかー．
形態の違いに目がいってしまうけど，
じつはその「行動様式・習性」がおおいにちがうのだ，
というご指摘もあったな．
-
真社会性昆虫にとって「コロニーそうじ」
はたいへん重要なお仕事
(私はこれと「(放浪性アリの) 分巣しやすさ」
みたいなことが関係あるんでは，
という妄想にとりつかれている)．
アブラムシは若い兵隊がそうじ担当で，
そのアタマに「ほうき」状構造があるとか．
-
そうじといえば，
コロニー内のアブラムシ死体は死後一日後ぐらいから
脂肪酸 (死亡酸とメモには誤記してしまった; リノール酸)
が「わいて」くるので，
それを識別してコロニー外に捨てられる．
シロアリだと (オレイン酸がわいて?) 死体は食われてしまうとか．
-
カースト認識も体表ワックスでできるものらしい．
-
コロニー内組成は個体間相互作用によって状況応答的に柔軟に変化する，
は誰にとっても受け入れやすいハナシになりつつあるようで
……
しかし個体の生理状態 (とその履歴)
とかは反映されないのかしらん?
たとえば母親の胎内にいるときや幼少時に飢餓状態を体験すると
エサ集め担当に分化しやすい，
とか．
うーむ
……
-
ところで，
この勉強会に参加したもひとつの理由は
……
「社会性昆虫を研究するヒトたち」
それ自体の観察だ．
まあ，
東さんのアヤしげな科研ぷろぢぇくとにまきこまれている以上，
このヒトたちを「ほほー」と言わせる計算 (計算?!)
を何かやってみせねばならぬわけで
……
といった「敵状」偵察ですね．
しかし誰も「計算」なんぞは必要としてないようにも見えるな
……
-
1830 ごろ終了．
いやー，
つかれました．
そしてまた熱帯林メイルとか
……
-
2015 研究室発．
2030 帰宅．
体重 72.2kg．
晩飯．
-
[今日の運動]
-
いかん，
積雪までの貴重な時間がうんどう不足の日々のうちに
……
-
[今日の食卓]
- 朝 (0830):
米麦 0.4 合．
チンゲンサイ・ハクサイ・豆腐のカレー．
- 昼 (1330):
研究室お茶部屋．
食パン．
- 晩 (2130):
米麦 0.7 合．
チンゲンサイ・ハクサイ・豆腐のカレー．
2006 年 11 月 09 日 (木)
-
0730 起床．
0745 自宅発．
曇．
そういや，
昨晩のラジオのニュースによると，
今朝は水星の太陽面通過が見られたそうな．
0800 研究室着．
朝飯．
コーヒー．
-
うだうだとぎょーむ日誌かきとか．
-
午前中は，
(多忙ぎょーむからいったん開放されてよーやく大学にくることできた)
志水さんと地衣類モデリング相談．
まあ，
だいたいこういうかんぢでよろしい，
ということで．
当方はアタマからこの階層ベイズモデルの詳細が消えぬうちに，
作文として残さねばならないな
……
-
熱帯林データどたばた問題，
関係者全員の総意として
当方に面倒をおしつける方向でハナシが収束しつつあるようだ
……
ということで 37.6 万個体 20 年ぶんのデータをながめて，
作戦を考える．
-
昼飯．
-
日没ごろまでこの莫大なデータから「使いものになる部分」
を選びだす Perl フィルターをくみあげる作業．
熱帯方面に問い合わせたり，
甲山さんにも採否の判断を投げたり．
-
生態学の莫大なデータというのは「まちがいのそうくつ」
みたいなものだが，
かくのごとく数百万アイテムにのぼるデータともなると，
もうどこもかしこもエラーだらけに見えてしまう．
簡単なやつはどんどんフィルタリングできるんだけど，
そのあげくに残存しているまちがいはなかなか手ごわい．
R
はまちがい発見能力に関しては Perl より格段にすぐれているので，
R にも手伝わせる
……
がこいつはその高性能とひきかえに動作がのろかったりすることもある．
-
けっきょく
-
2005 年観測時の新規加入個体は対象としない
(新規入力データは不完全かつまちがいだらけなので)
-
2000 年観測値が欠測な個体は使わない
(甲山さんのよくわからぬ趣味)
-
観測はされているけれど，
2005 年時点のコメントで D (死亡) または B (破壊)
とあるものは除外
(2000-2005 年のあいだの D/B 個体数は 2.3 万ぐらい)
-
バットレスのせいで 2005 年どうまわりが測定できてない
ちょー巨大個体はいちおう採用
(2000 年観測値で一部を代用 …… これまた甲山さんの趣味)
-
ヘンな座標の個体もとうぜん除外
と選抜すると 25.7 万個体まで削減された．
-
お，
マレイシアからも「絶対必要です」攻撃が．
-
さらに細かい選抜がつづく．
この段階で残存・潜伏しているのはかなり「独創的な」エラーだったり．
-
とりあえず R で作業できるところまでスクリーニングが終了．
まだまだ長びきそうなので晩飯調達してきてお茶部屋で晩飯．
送粉ねっとわーくとかを説明してる文献勉強してる院生たちと
connectivity なる index を検討
……
なんかこれってアヤしげな割り算値だという気がする．
まあ，
グラフ上の統計量は組合せ論的というか数えあげ値の割り算値になりがちだけど．
-
タイムリミットつき
国際救助隊
ぎょーむ，
まだまだつづく．
このあとはいよいよサンダーバード 2 号で現場に輸送された救助メカならぬ
以前つくった自動作表・自動作図の R プログラムたちを投入するばかり，
なんだが
……
前回の (2000 年データと詐称していたけどじつは) 1995 年データと
今回の 2005 年データはあれこれとかなり異なる
(これは前回のデータが人手をわたるうちにこねくりまわされたため)．
ということで，
救助メカたちを活躍させるための状況，
「書式」みたいなものをととのえる作業がいろいろあるわけで．
「絶対必要です」仕様変更もあるしね．
-
なンかまたあとでやらされそうなので，
その手順を記録．
-
SPRITdata.R
で 50ha データを樹種別・セクター別に分割
(これは元データが変更されたときに一回やればよい)
-
select_all.R
で 50 樹種 1000 個体を 9 set 選抜．
-
make ds で
generate_datasheet.R
が動いて LaTeX ファイルを出力，
その後に LaTeX が召喚されて 50
ペイジの野外調査シートの PDF ファイルを生成
-
make maps
(特にこちらは R のフォントまわりの面倒を回避するため，
R --vanilla モードで動かす
make を使う)
で generate_map.R
が動いて標的 1000 個体ぶんの近傍 5m 四方の局所地図
84 ペイジぶんを R から直接 PDF 出力;
さらに
map50ha.R
で 50ha 全図も出力
-
作業終了 2315．
アップロードは明朝でいいだろう．
(これは毎日やってることだが)
差分ばっくあっぷをとり
2325 研究室発．
雨．
2340 帰宅．
ばてた．
-
[今日の運動]
-
[今日の食卓]
- 朝 (0820):
研究室お茶部屋．
食パン．
- 昼 (1310):
研究室お茶部屋．
食パン．
- 晩 (1930):
セイコーマートのサンマ鮨．
380 円．
2006 年 11 月 10 日 (金)
-
0745 起床．
0755 自宅発．
曇．
0810 研究室着．
朝飯．
コーヒー．
-
とりあえず昨晩生成したマレイシア調査用 PDF ファイルのうち，
必要最小限のものだけ選抜して圧縮．
それでも 1.4MB あるか．
で，
ダウンロードできるところにおく．
そうか時差が一時間あったな．
-
ぎょーむ日誌とか昨日から放置してるメイル返信かきあれこれとか．
また査読しやがれメイルとかもあるし．
-
R のこれまたよくある質問．
グラフの軸で単位とかかきたいときのベキ乗 (上つき)
とかはどう書けばよいか．
ぎょーむ日誌でも
何度か紹介してるけど，
たとえば
expression("density (g cm" ^ -3 * ")")
とすれば，
となる．
expression()
の
() 内に
「数字だけでなく文字列をも項とする疑似数式みたいなものを書く，
二項演算子でつなぎつつ」
と考えるところがポイントである．
すると * は「かける」をあらわし
expression() はかけ算記号が省略した表記
("x" * "y" はかけ算記号を省略して "xy" と表記する，
といったルールになっている)
で graphic device 上に描くという対応関係が理解できるだろう．
-
また院生のマシンに
R のインストール
& Tinn-R
のインストール．
-
アップロードしたあれこれに関して，
マレイシアからはやりなおし命令もなく「問題なし」との連絡が．
IR (いんたーなしょなるれすきゅー)
ぎょーむ，
無事に終了，
か?
-
12/7 (木) - 12/10 (日)
東京出張
の宿の手配．
昨年 12 月
の R ゆーざー会のときは，
羽田から至近の蒲田に泊まった，
というのも南麻布
には土曜日の朝に移動すればよかったので．
しかし今回は金曜日朝の通勤通学時間帯にも移動せんといかんわけで
……
蒲田から都心への経路はスシづめ状態だろう．
ということで軟弱にも統数研から至近距離に宿をとる
……
一泊 9450 円!
うーむ
……
蒲田あたりだと 6900 円ぐらいなンだが．
航空券も予約．
-
竹中さんから
5th International Workshop on Functional-Structural Plant Models
(FSPM, 2007 年 11 月 4-9 開催)
の案内いただいたわけだが
……
いやはや，
要旨提出は来年の 3/19 か．
ぱいぷ樹木
の改良改善にかまけているヒマあるんだったら，
アリ科研費しがらみ研究と未提出宿題あれこれやりなさいよ，
という「借金」状況だよなあ
……
ああ，
にゅーじーらんど行ってみたかった．
-
だったらアリ関係のメイルでも書けよ，
ということになるが
……
とりあえずすでに午後だし，
気力もつきたので昼飯．
-
MacOS X のヒトたちも
Darwine
とかインストールすれば WinBUGS が動かせる (??)
だろうから
R2WinBUGS
とか使えるようになるのではなかろうか
……?
特に Intel Mac 上ではかなり高速に動作するような気がするんだけど．
-
といったことを考えつつメイルなど書いたり．
しかしアリ関連メイルのほうは書く気合いが足りない．
-
なンかいろいろ雑用とか．
-
そして小林さんの新アイデアでアカマツ (野外データのほう)
解析やりなおし
……
うーむ，
やはり「個体差」とか考慮せんといかんのか?
[光合成速度 vs 葉重量]
まあ，
最尤推定値は「まがる」わけですが．
AIC で比較すると，
単なる直線回帰が勝つんだよね．
そしてわれながらバカげていることやってみたんだが
……
「一個体から一データしか取れてないけど『個体差』推定しろ」
とかいう命令だしても，
R
の
lmer() や
nlme()
はなかなかカシこくて
(というか当然ながらというか)
「そんなのできるか」
と却下されてしまった．
これまた順当なことにというか，
処理ごとにセンをあてはめる，
というのもパラメーター使いすぎで AIC よくない．
-
1900 研究室発．
1920 帰宅．
体重 72.4kg．
晩飯．
すぐ寝てしまった．
-
[今日の運動]
-
[今日の食卓]
- 朝 (0820):
研究室お茶部屋．
クロワッサン．
- 昼 (1320):
研究室お茶部屋．
食パン．
- 晩 (2040):
米麦 0.7 合．
ダイコン・豆腐・ワカメ・煮干の味噌汁．
