ぎょーむ日誌 2007-12-(21-31)
2007 年 12 月 21 日 (金)
-
0730 起床．
ねむい．
コーヒー．
朝飯．
0855 自宅発．
晴．
0910 研究室着．
-
神山さん修論の光計算のみなおし
……
うーむと思っていると，
神山さんから階層ベイズモデルの補足説明を，
ということでまずはそのあたりを「脱」割算値計算にそって
最初から説明しなおしてみる．
で，
そのあと光計算に関してうーむとナヤんでみる．
うーむ．
-
散乱された光は必ず「下」にいく，
といったいくつかのわかりにくい仮定が潜在しているのしらん?
その他を納得することにして透過・散乱・吸収の過程の部分だけ，
神山さんとひととおり整理してみる．
その後，
昼飯に．
-
また光計算みなおし作業．
やはりあちこちにうーむな部分が
……
-
途中で一時間ほど石狩浜海浜植物作図問題，
本日は「あの」ハマエンドウ作図ふたたび
(すでに今年の 2 月初めにハマエンドウ個体ごとの葉食害作図で
苦闘
ずみ)，
である．
今回は個体ごとの繁殖とシュート成長の同時表示で，
こいつは花序が複数あったりなかったり，
分枝があったりなかったりでめんどうきわまりないもの．
さらにたいへん申しわけなさそうに
すごく
難しい条件を追加したりするのが
稲葉さんの特技だったりなかったり
……
-
なんとかなったので，
また神山さんと光計算まにあっく談話．
この
Beer-Lambert
なるちょっと独特なる空想的な
世界では吸光係数とか呼ばれることの多いパラメーターに
およそありとあらゆる面倒という面倒をことごとく押しつけるので，
その他の部分はもうちょっとすっきりするはず，
というような．
そして，
そもそもゑくせるで計算できる範囲からすごく逸脱してるような
……
-
1800 前に神山さんが仙台に帰られる．
で，
1800 からまた別のデータ解析こんさるが入るんだが
……
その前に PC が起動しないので修理してほしい，
という院生ハードウェアこんさる．
これは「はんどぱわー」で 10 秒ほどで修理できた．
私にこんな秘められた才能があったとは．
-
1800 すぎから生態遺伝の岡田君のデータ解析こんさる．
「折れ線」回帰か
……
2003-10-29
のコードを参考にしようとしたんだけど，
ほとんど作りなおすことに．
なぜかしら「折れ線」のほうが格段に AIC とか低いなぁ
……
と不思議に思っていたんだけど，
原因はこのデータがいわゆる「等分散」ではないことにある，
と気づいた．
データのばらつきに gradient がある場合には，
(線分ごとに自由に分散を選べるタイプの)
「折れ線」回帰が対数尤度向上に関してはすごく有利，
ということだ．
こういう統計モデリングがいいのかわるいのかよくわからんけど．
-
意外と苦闘してしまって 2 時間以上をついやしてしまった．
その後はまた稲葉さんと石狩浜作図問題．
ハマエンドウ 120 個体の葉数変化，
library(lattice)
で個体ごと一覧を生成．
-
2100 研究室発．
買いもの．
2120 帰宅．
ぐったり．
晩飯の準備．
晩飯．
-
[今日の運動]
-
[今日の食卓]
- 朝 (0750):
食パン．
- 昼 (1250):
きゃら亭．
今日の A 定食はサフランライスのシーフードピラフ．
- 晩 (2310):
米麦 0.8 合．
シイタケ・ワカメ・豆腐の味噌汁．
キムチ．
2007 年 12 月 22 日 (土)
2007 年 12 月 23 日 (日)
-
0910 起床．
朝飯．
コーヒー．
洗濯．
怠業．
-
うーむ，
なぜか東部ニューギニア戦史
(Wikipedia:
日本語版,
敵国語版;
公刊戦史かしらん?
Victory in Papua
← 地図がくわしい)
をだらだらと調べて昼すぎになってしまった．
人間は
こんな
ところでも戦争できるんだなぁ
……
というか，
そもそもこんなところに以前から原住のヒトたちが集落を作ったり標高 2000 m
より高い山越えの道ができてたりするんだなぁ，
というほうがすごいか．
-
1315 自宅発．
曇．
1330 研究室着．
昨日のぎょーむ日誌などうめてみる．
アップロードしたあとに気づいたんだけど，
昨日はヘンな査読依頼なんかもあって仕事がススまなかったんだ．
えーい，
これどうしてくれよう．
しばらく放置．
-
田辺さんと苫小牧 & 屋久島 33 樹種の実生の統計モデリング改良の相談．
えーと，
種子重の効果と
「 pioneer かどうか」効果をみる，
と．
これは説明変数にする方式と，
分散共分散行列で表現する方法の両方があるんだけど，
どうしたもんかな
……
「明るい場所」
うんぬんは種差と区別がつかない，
低木・針葉樹は種数が少ないのでそういうよけーなカテゴリーわけは不用．
-
しかし「苫小牧か屋久島か」なる「樹種の上のカテゴリー」
なるものを導入させられそうなかんじで
……
これもあとまわしだ．
M2 修論 5 つのうち，
樹木種数多数モノつまり甲山さん的研究が 3 つあるんだけど，
この甲山さんの意図がなかなか読めなくてモデリング方針が決めにくい．
-
とりあえず宮田さんの苫小牧樹木てっぺんシュートのモデリングの中で
つかんだ「巧妙なる個体差導入による deterministic node
連鎖切断による MCMC 計算高速化」
によってコードを書き換える
……
ついでに抜本的なコード整理をしておくことで，
上のようなごちゃごちゃした悪しきモデル拡張にそなえる．
-
いろいろと試行錯誤．
とりあえず苫小牧樹木だけで試験運転
……
2000 MCMC step に 450 秒ぐらい．
以前に比べて数倍速くなり，
収束もだんぜんよくなった．
しかし現時点では，
まだ上のごちゃごちゃ要因もいれてないし，
「種差」パラメーターの事前分布は独立した正規分布のままである
……
これはあとで多変量正規分布に置換せねばならない．
-
2125 研究室発．
2140 帰宅．
晩飯．
-
[今日の運動]
-
[今日の食卓]
- 朝 (0930):
食パン．
- 昼 (1330):
研究室お茶部屋．
米麦 0.6 合．
海藻スープ．
- 晩 (2230):
米麦 0.8 合．
ジャガイモ・ニンジン・タマネギ・豆腐の味噌汁．
キュウリ．
2007 年 12 月 24 日 (月)
-
0810 起床．
また電気つけっぱなしで寝てしまった．
なんかばててるなぁ．
朝飯．
コーヒー．
なんだか気分がいまいちよろしくなくてうだうだする．
1140 自宅発．
雪．
1155 研究室着．
-
昨日の統計モデリング改良
……
というか複雑化のつづき．
とりあえず種子重の効果をいれてみる．
やはり収束が遅くなるな．
ということで burn-in 期間を延伸する必要がある．
とうぜんながら，
よけーな変数とパラメーターがあれこれ増えたので，
計算そのものも遅くなる．
-
計算まち時間に昼飯．
-
個体ごとの種子重なる未観測値があまりばらつかぬよう調整したり，
初期値を丁寧に設定したり，
非線形な部分をへらしたり，
library(lattice)
図で収束の悪いことに関係ありそうな説明変数を見直したり
(下の図は
xyplot(log(Hs) ~ log(Wsl) | Sp, data = d, layout = c(6, 2)))
して統計モデリングの改善，
こういう地味な修繕の積み重ねによって MCMC 計算は高速化され収束も改善される，
と
……
-
そして独立な事前分布から生成される事後分布の
pairs()
図とはかくのごときものだ．
うーむ，
けっこうな労力をついやして多変量正規分布な事前分布を導入すべきかどうか
……
とりあえず，
とうぶんあとまわしにすることにしよう．
-
さて次に「耐陰性」なるかなり主観的な変数をどうくみこむか，
だけど
……
このあたりベイズ的な工夫をいろいろやってみたんだけど，
どれもいまいちぴんとこない．
というか (形式的には) random effects 的に入ってる「種差」
なるパラメーターとほとんど同じ挙動，
といいますか．
-
けっきょく，
このあたりは説明変数として乱暴に投入．
まあ，
そもそも乱暴なグループわけ法なんだからこれでよいのかもしれない．
おそらく，
樹種の「上」に存在するグループわけとしてはほとんど機能しないんじゃないかしらん?
-
こういう乱暴モデリングだと
2000 MCMC step 420 秒で
……
うーむ，
収束はいまいち．
そしてこのいいかげんな説明変数がいくつかの資源分配・形態形成を説明している?
-
本日はどうやらお茶部屋における大学院生の活性がいつもよりさらに高まる
特異日
のようだ
……
とういことで，
嵐にまきこまれぬよう早々に退避することに．
1810 研究室発．
ひさしぶりに札幌駅周辺の本屋あちこちをふらふらとハシゴしてみる．
2000 帰宅．
晩飯の準備．
晩飯．
-
[今日の運動]
-
[今日の食卓]
- 朝 (0900):
米麦 0.6 合．
納豆．
- 昼 (1330):
研究室お茶部屋．
米麦 0.6 合．
- 晩 (2030):
米麦 0.7 合．
ジャガイモ・ニンジン・タマネギ・豆腐の味噌汁．
2007 年 12 月 25 日 (火)
-
0820 起床．
朝飯．
コーヒー．
0940 自宅発．
曇ちょっと雪．
積雪量，
今年はまだまだ，
というかんじで．
まあちょっと大雪になればたちまちにして歩道と車道のあいだに
巨大な除雪による雪壁が形成されるわけだが．
0955 研究室着．
-
苫小牧 & 屋久島実生データのうち，
とりあえず苫小牧のほうのモデリングを改善する試行錯誤のつづき．
-
サロベツ種子発芽実験データ解析でもちいたやや変則的な Bayesian ZIP モデルの
説明作文かきの手つだい．
院生たちは確率なる概念を道具として使えるけれど，
確率変数 (random variable) についてはもう少しあいまいに理解しているようで．
-
実生データの統計モデリング，
「耐陰性」なるアヤしげな変数をベイズモデル中の確率変数
(random variable) としてあつかうための試行錯誤．
屋久島データとは異なり，
苫小牧データには「耐陰性なし」「耐陰性あるのかどうかよくわからん」
の 2 水準しかないせいか，
後者はことごとく「耐陰性あり」に assign されて動かなくなった．
そのくせ全体の収束は遅い
……
全体というか，
なぜか重量の測定誤差パラメーターの収束がゆっくりになっていて，
全 deviance もそれにつられてゆっくり変化．
-
そしてこの「耐陰性」なるじつにアヤしげな説明変数が，
苫小牧 13 樹種を 5:8 に分割してるんだけど，
なぜかしらけっこう実生のカタチを説明する「グループわけ」になっている，
というおどろくべき推定計算結果が
……
うーむ．
-
現時点では 4000 MCMC step で 990 秒ぐらい．
収束改善のため，
もっと計算時間を増やすことは可能なんだが
……
-
これまた面倒なことに，
苫小牧 & 屋久島データを同時にあつかう，
といった統計モデリングが必要かもしれない，
ということでそれにあわせた改造を．
まあ，
その準備はしてたんで問題なし．
試験運転しようとすると，
屋久島のナゾの樹種コード
Sp
(4 個体だけ)
の先験的「耐陰性の値」がないことに気づいた
……
えーい，
この樹種はいったん削除だ．
-
試験運転．
昼飯．
-
なんというべきか，
くだんの「耐陰性」なる確率論的なグループわけ，
苫小牧と屋久島ではまったく別の方向に効果があるようなないような
……
「これ以上ややこしくするのはイヤだなぁ」
と思いつつも，
耐陰性と屋久島の交互作用項じみたものを追加してふたたび試験運転．
やれやれ．
-
えーい，
なんだかおかしいと思ったら，
この先験的なる耐陰性なる labeling の意味を逆にとっていた．
ラベル 1 は「耐陰性が低い」ではなく「高い」のか!
いやはや．
ミズナラが耐陰性ひくくて，
バリバリノキが耐陰性たかいのはいくらなんでもおかしいもんなあ．
樹種名が
Qm とか La
とか略されていたもんで
……
(このハナシはまたあとでひっくりかえる)
-
で，
逆にすると「何で?」とゆーぐらい収束が改善された
……
あいかわらず測定誤差のばらつきパラメーターは遅いけど．
で，
逆にしてもやっぱり苫小牧と屋久島では「耐陰性」なるラベルの働きかたが
ちがうみたいだなぁ
……
一例をあげてみると，
苫小牧では「耐陰」な連中は地上部重視だけど，
屋久島では地下部重視とか．
ホントかな?
-
まあ，
この観測データはややこしい点も多々あるから
……
-
このあたり，
またモデリング修正を余儀なくされる可能性大だけど，
もうちょい MCMC サンプリングを長くしてみるか
……
4000 MCMC step に 1890 秒．
収束はいまいちかなぁ
……
-
といろいろ検討してると，
データを入力した田辺さんがきて，
tolerancy 列でミズナラやハリギリが「耐陰性ナシ」になってるのは入力まちがい，
といった修正を
……
うーむ．
ということで，
また MCMC 計算やりなおし．
-
2000 MCMC step 953 秒
……
で，
もとデータを修正するとこれまた「何で?何で?」というぐらい収束が改善してしまった．
うん，
まあ改善されたんなら何でもいいや．
(事後分布表，
BUGS code)．
さっきとは逆で，
苫小牧では「耐陰」な連中は地下部重視だけど，
屋久島では地上部重視になってしまった．
-
さて，
これで苫小牧・屋久島データの統計モデリングがまとめてできていればよいのだが
……
作図プログラムはかなり改善せんといかんな．
とりあえず，
こんなかんじとか?
-
Yak + Tol + Yak:Tol
なる交互作用のいれかたはまちがいではないけど解釈がめんどうになるところがあるな
……
ということでそのあたり修正して計算やりなおしを命じて本日は撤退．
まあ，
苫小牧 + 屋久島データ統合して 15 分ちょいで推定計算が終了するようになったのは，
おおいなる進歩だ．
1930 研究室発．
買いもの．
1950 帰宅．
晩飯の準備．
晩飯．
-
[今日の運動]
-
腹筋運動 30 ×
3 回．
腕立ふせ 10 ×
3 回．
スクワット 100 回．
-
[今日の食卓]
- 朝 (0850):
米麦 0.7 合．
ネギ卵炒飯．
- 昼 (1300):
研究室計算機前．
食パン．
- 晩 (2240):
米麦 0.8 合．
ワカメ・豆腐の味噌汁．
ハクサイ・ニンジン・タマネギ・ネギ・ピーマン・ブナシメジ・イカの炒めもの．
2007 年 12 月 26 日 (水)
-
0820 起床．
さいきん「電灯つけっぱなし寝」が増えてしまっている
……
朝飯．
コーヒー．
0940 自宅発．
晴．
いつのまにか歩道が雪で被覆されている．
0955 研究室着．
-
苫小牧 + 屋久島統合モデリングのつづき．
意味のない中央化をやってたことによーやく気づいたので削除
……
しかしあとになって，
まあこのへんはどう「ずらし」ても結果は同じだし，
中央化してたほうがゐんばぐすの収束は良いのでそのままに．
昨日わずらわされていた「耐陰性」なる確率論的なグループわけは
初期値のおきかたにけっこう依存することがわかってきた．
このあたり修正しつつ，
再計算をくりかえす．
-
いろいろとややこしい推定結果がでてくるので，
いまさらながら観測データ作図を工夫してみる．
うーむ．
なンというか，
べいづな方法は強力なのでかかるノイズだらけの中から
何か傾向みたいなものをあらわす事後分布をすくいだしてくるわけだが
……
2007 年 12 月 27 日 (木)
-
0820 起床．
朝飯．
いろいろとあって出発がおくれる．
1110 自宅発．
曇．
1125 研究室着．
-
苫小牧 + 屋久島 29 樹種実生の統合モデル，
4000 MCMC step に 2000 秒を費やして
……
まあ，
こんなもんでしょう，
という収束
(事後分布表，
BUGS code)．
-
作図プログラムを工夫．
まずは
「切片」的なパラメーター
……
というか「種間共通 + 種差，
ただし種子重量の効果と『耐陰性』ラベリングの効果をのぞく」
といったところか．
-
「種差」事後分布の
pairs() 図．
とうぜんながらというか，
いっそう「丸く」なった．
-
で，
サイズ・種子重そして「耐陰性」ラベルつけの効果．
-
とりあえず昼飯．
-
統計モデルの挙動を一覧する R 作図コード書き
……
データ構造のややこしさ
+
library(lattice) のわかりにくさ，
で苦闘．
-
夕方までぢたばたしてよーやくできた．
6 種類の数量的関係に関して，
たとえばこのように全 29 樹種に関して
……
-
以上の結果を，
田辺さんに納品
……
これで M2 下うけ 5 件中 4 件はだいたい終了．
脱力．
-
5 件目はしばらく着手しなくてよさそうだ
……
2040 研究室発．
米 10 kg 買う．
3020 円．
2100 帰宅．
体重 70.2 kg．
晩飯．
-
脱力のあまり，
ちょっと風邪ぎみ?
というほどでもないな．
よくウガイすることにしよう．
-
[今日の運動]
-
[今日の食卓]
- 朝 (1030):
米麦 0.7 合．
ネギ・豆腐の味噌汁．
カスベ
電磁波酒蒸し．
キャベツ・ニンジン・キュウリ・タマネギ・ハムのサラダ．
- 昼 (1400):
研究室お茶部屋．
米麦 0.6 合．
春雨スープ．
- 晩 (2230):
米麦 0.8 合．
ハクサイ・豆腐の味噌汁．
キャベツ・ニンジン・キュウリ・タマネギ・ハムのサラダ．
2007 年 12 月 28 日 (金)
-
0840 起床．
体力的には回復，
かしらん．
朝飯．
コーヒー．
なぜか洗濯．
1010 自宅発．
曇．
1025 研究室着．
-
M2 したうけ奔流がしばらく止まることになった
……
のでいろいろとこの間に放置してたあれこれを片づけねばならないんだけど
……
ちょっと気力ぎれ．
-
昨日の計算など確認してると，
(ここしばらくの面倒な統計モデリングそもそもの根源，
とちまたでウワサされている)
甲山さんがうろうろと．
「多樹種」なあれこれに関する議論
……
うーむ，
オトシどころについてはやはりあまり考えては
……
-
昼飯調達のため北大生協往還の旅．
積雪が少ない．
書籍部もうろうろ．
あまりにもぼけーとしていたので，
ノドあめを買うのを忘れてた．
まあ，
そんなに深刻にヤラれているわけではないが．
研究室にもどって昼飯．
-
午後もぐったりぎみ．
お，
Vine Linux
の 4.2 がリリースされた．
ということで，
さっそく
/etc/apt/sources.list
を 4.2 対応にして
sudo apt-get update && sudo apt-get dist-upgrade
……
glibc の Vine patch version が変わったぐらい?
あまり package 更新がない，
というか．
-
で終わるのはもったいないので放置雑用にとりくんでみる．
-
生態学会さーばー雑用．
これは一瞬でおわった．
さっさとやればよかった．
-
メイルかきあれこれ．
2/14 (木)・2/15 (金)
はつくばの森林総研で「ばれんたいん」統計学講義
(出張期間は 2/13-2/15)．
といっても内容は
10-11 月の授業
と同じにするつもりだけど．
-
考えてみれば森林総研には
伊東宏樹さん
がおられるんだから
(ただしつくばではなく多摩)，
私なんかが引きうける必要なかったような
……
今度からこういう依頼のときは先方の状況をよくみることにしよう．
-
よくわからぬ
「GCOE 申請したいんだよー」
関連ファイル．
前回敗退したのでまた出すんだ!
とかゆーことらしいんだけど
……
何やら私にもぎょーせきファイルの提出が要求されている．
なぜ無関係なはずの私がかくのごとく巻きこまれてるかと言えば，
(1) 環境科学研究員の教員であり
(2) 2004-2007 年の期間に院生との共著論文があるから
(浦口さんのカエデ論文)．
-
GCOE では「教育」とやらがじゅーしされてるとのことで，
大学院生と教員が一致協力して論文をだしました，
といったウルわしいハナシの個数が何かを represent
しているのだッということになってるらしい
……
嗚呼，
修論の論文化といえば岩倉さんアリ論文とか牛原さん屋久島葉寿命論文とか
……
このあたり，
「真の」指導教員のかたがたにある種の問題があるよーな気も．
-
で，
送られてきたのはわりと兇悪なゑくせるファイルであり
……
しかも Shift-JIS 文字コードを使ったファイル名を変更すんな，
とゆーナイスすぎる拘束条件つき
(とゆーことで放置されてた)．
ふと思いついて，
Linux 上の VMware のゐんどーづ上の web browser で Gmail
を使ってくだんのファイルを download，
ゐんどーづ版の
OpenOffice.org
を使ったらこの罫線まみれの呪われきったファイルが編集できた．
そのまままた Gmail 経由で提出．
-
しかし北大もいいかげんに事務処理をもっとデータベイスと連係しろよ，
と思うわけだが
……
まあ，
とうぶんムリなんでしょうなぁ．
紙書類ベイスの情報処理時代と比較して，
紙がゑくせるに置き換わっただけだ．
ぎょーせきデータとか北大データベイスに入ってるんで，
そのあたりをいくらでも工夫できるはずなんだけどねえ．
やれやれ．
理想的には大学事務処理用の
「身軽な」プログラミング担当な技官さんたちがいればよいのだが．
時代錯誤ぎみの大型計算機「すぱこん」のおもりを自慢してないで，
もっとネットとかデータベイスプログラミングとかちゃんとやればいいのにねえ．
いまやより重要なのは，
ハードウェアよりソフトウェア，
ソフトウェアよりデータなんだってば．
-
VMware つけたついでに
……
先日，
ef-prime
の鈴木さんからメイルいただいたので，
R Analytic Flow
のプレヴュー版を試してみることに
(現時点では Windows 版しかないので)．
Java まわりのインストールが必要．
下は Linux 上の VMware 上のゐんどーづ上で R Analytic Flow
が動作してる様子．
-
予想してたより R が「前面にでてる」ソフトウェアになっている
……
おそらく (現 version では)
R の知識がないヒトが
「解析の流れ」のノードをぽちぽちと追加して解析 tree
を構築ことは不可能ではないかな．
むしろ使いかたとしては
(ここしばらくの私のぎょーむに照らして言うと)
「下うけ」解析の結果を「納品」するときにこれで
データの流れを表現して，
ノードをぽちぽちと選択して「実行」すれば
「ここがこうなってあーなってこういう結果が得られるんですよ」
てなことを示せる
(あるノードを選択して実行するとそこから「上流」の行程が
R の中で実行される)，
ってことかな．
-
データ解析とあまり縁のないお客さん相手にはなかなか有効かも
……
というか，
よけーなグラフィカルインターフェイス (GUI) がついてない R
だからこそこういう Java ソフトウェアが作れたりするのかもね．
-
私なんかも客商売だけど，
主たる相手が大学院生なので問答無用で文字ばかりの
CUI
を平然と押しつけさせておりますが
……
そして「計算機といえば GUI」世代なヒトたちも実は
CUI なんかを平然と使いこなせるようになるわけですが．
-
それはともかくデータ解析会社 ef-prime
はますます繁盛しておられるようでめでたいことです．
-
粕谷さん科研費を使わせていただいて，
出張するための連絡メイルを九大の生物事務室に．
-
さーて，
次なる M2 下うけがいつ発生するのかは神のみぞ知る，
という状況のもとで，
残された雑っぽいことは先日きた不幸なる査読依頼を
うまく断ること
……
なぜ断りたいかといえば投稿者があまりマトモではなさそうなヒトで，
査読者・編集者が何度も指摘してるのに，
かなりおかしい手法に固執してぜんぜん修正する気配がないため．
-
あとアヤしげなデータ解析メイルが放置しているな．
それから雑っぽいことではないけど，
志水さんの地衣類論文の投稿前最終みなおしと，
シウリザクラ論文の再投稿のためのみなおしがある，
と．
-
それからアリ関連だと
……
アリ論文をさっさと書きたいんだけど，
その前に生態学会のアリポスターの要旨を準備して皆さんに
連絡しないといけないな (提出〆切は 1/8)．
-
そうそう，
生態学会福岡大会，
プログラムがきまったようで．
えーと，
ちゃんと定刻どおりに出現すべきは
……
てなところかな．
前半に集中．
3/16 (日) 15:45 の ECOAS 企画集会も出とかないといけないんだろうなぁ．
-
わんぽいんと M2 リリーフなんかもあったりして，
1950 研究室発．
2005 帰宅．
晩飯の準備．
晩飯．
-
[今日の運動]
-
[今日の食卓]
- 朝 (0920):
米麦 0.6 合．
ハクサイ・豆腐の味噌汁．
- 昼 (1330):
研究室お茶部屋．
北大生協にぎり飯．
- 晩 (2120):
米麦 0.7 合．
ハクサイ・油あげの味噌汁．
ピーマン・卵の炒めもの．
2007 年 12 月 29 日 (土)
-
0940 起床．
すごーくよく寝た
……
ノドの調子なんかもよくなっている．
朝飯．
コーヒー．
怠業．
-
昨日の Vine Linux update で
TrueType-sazanami
のたぐいがインストールされたので，
とりあえず sudo apt-get remove TrueType-base．
かわりに ver.2.5 ぐらいの古い Vine
置場 (FTP site) から
TrueType-base
&
TrueType-kochi
をとってきて rpm -Uvh
……
ここで RPM データベイスをごまかすために
sudo rpm -e --justdb TrueType-base TrueType-kochi
する必要ある．
で，
このあと下のごとく
$HOME/.Xresources
を変えてみた．
Xft*antialias: True ! ここを True にすればアンチエイリアスがオン。
Xft*dpi: 96 ! GNOME ではこの値にしてるようなんで、合わせといた方がいいかも。
! こっから下は、アンチエイリアスをオンにした場合の設定。
Xft*hinting: True ! 多分これは True の方がキレイ。
Xft*hintstyle: hintfull ! 他に hintnone, hintslight, hintmedium というのがある模様。
Xft*rgba: rgb ! サブピクセルレンダリングの設定。詳しくは GNOME のコンパネで。
Kinput2*background: #aaaaaa
Kinput2*OverTheSpotConversion.modeBorderForeground: false
Kinput2*OverTheSpotConversion.modeLocation: bottomleft
Kinput2*Font: -alias-fixed-medium-r-normal--16-*
Kinput2*KanjiFont: -alias-fixed-medium-r-normal-*-16-*-*-*-*-*-jisx0213.2000-*
Kinput2*KanaFont: -alias-fixed-medium-r-normal-*-16-*-*-*-*-*-jisx0213.2000-*
KTerm*savedLines: 2048
KTerm*background: #001000
KTerm*foreground: #40e040
KTerm*cursorColor: #e0a040
KTerm*scrollBar: false
KTerm*VT100*fontList: -alias-fixed-medium-r-normal--16-*
KTerm*VT100*boldFontList: -alias-fixed-bold-r-normal--16-*
みればわかるとーり，
私はいまだに古くさい kinput2 と kterm
に依存している．
anthy + SCIM
を使おうとときどき努力してみるんだけど
何だかいろいろと気に入らないので kinput2 + Canna
にもどる．
-
Terminal に関しては
kterm
がもっともすっきりしている．
しかし最近は部分的に gnome-terminal
を使うようにしている．
今日しらべていてわかったんだけど，
こいつの文字表示の設定にはコツがあり，
(私はふだんは GNOME 環境なんぞは使ってはいないんだけど)
一時的に gnome-panel
でゐんどーづ的な状態にして，
上のめぬーの
[デスクトップ] → [設定]→ [フォント]
で
「フォントのレンダリング」
が default では「最適なシェープ」になっているので
(これはかなり anti ailiasing がかかっている)，
これを上の図のごとく「最適なコントラスト」
にすると
(私にとっては)
「見やすい」状態になった．
-
今月のガス代 7440 円
……
これは寒くない月の倍以上の請求．
平年より高いように思う．
プロパンガス値上げの影響か?
-
1220 自宅発．
雨．
なんでこの時期に雨なんだよー．
まあ毎年 180 億円も除雪費用に使ってる札幌市的にはありがたいのかも．
しかしなンと言いますか
……
でこぼこになった雪面によってあちこち水たまりが形成されてて歩きにくい．
というのも私が 2 年ほど前に Homac で買った安もの偽トレッキングシューズふう
の靴は安ものにふさわしく (つい先日なんだけど)
側面に破孔というか断裂が生じてしまって，
ですね
……
-
1240 北大構内クラーク会館地階の
理髪店
(なぜかリンクぎれ)
着．
1300 散髪開始．
前回のさんぱつは
10/10
だったのか
……
どたばたしてたとはいえ 2.5 ヵ月以上も放置してた．
やや短く切ってもらった．
すっきりさっぱり．
1400 研究室着．
昼飯．
-
Vine update の余波つづく．
私の X 端末機として使われてる
lex
(闇ネットの DNS server という大任も果たしているけど)
の NVIDIA driver の挙動がいまいち
……
そしてできあいカーネルでは driver のコンパイルもできん，
……
とゆーことでまた kernel コンパイル．
いやはや．
あ，
どうせなら新しいカーネルソースとってくればよかった．
-
以上の悪戦苦闘はじつにばかげた失敗に起因するもの，
と判明した．
lex
の video card は古井 GeForce 4 なのに GeForce 8
用の driver をコンパイルしてた．
互換性はまったくない．
正しくは Nvidia の driver download page に行って，
``Legacy'' なるカテゴリーの中から GeForce 4 用ソースをダウンロードせよ，
と．
この最新版は
NVIDIA-Linux-x86-96.43.01-pkg1.run
だった．
-
一件落着したので．
1900 研究室発．
1920 帰宅．
晩飯の準備．
晩飯．
いろいろとそうじというか洗濯というか
……
-
[今日の運動]
-
[今日の食卓]
- 朝 (1030):
米麦 0.6 合．
ハクサイ・油あげの味噌汁．
- 昼 (1410):
研究室お茶部屋．
北大正門前の
ぱん吉
で買ったフランスパン．
リンゴ．
- 晩 (2110):
米麦 0.7 合．
ダイコン・ニンジン・ジャガイモの味噌汁．
ゆでホウレンソウ．
ダイコンの酢のもの．
コマイの干物．
2007 年 12 月 30 日 (日)
[確率密度関数の憶測]
たとえばいま平均ゼロで標準偏差 1
の正規分布 (破線) から 100 個のデータが得られた
としよう．
このデータは
|
で示されている．
この
|
たちが「正規分布から得られた」
ということがわかっていればもともとの確率分布を推定するのは
まったく簡単．
しかし「どんな分布かわからん!」
という状況ではどうしたらよいか?
そういうときでも
R の
density()
を使って
実線で示しているような憶測確率密度関数が
(推定というより) 憶測できたりする．
-
えーと上の憶測計算の手順としては，
たとえばこういうかんぢで．
n <- 100 # データ数 n = 100 (個)
data.x <- rnorm(n) # N(0, 1) の乱数
dxy <- density(data.x) # density() を使った推定
plot(dxy$x, dnorm(dxy$x), type = "l", lty = 2) # もとの確率分布の表示
points(data.x, rep(0, n), pch = "|") # data.x の "|" 表示
lines(dxy$x, dxy$y) # 憶測された確率密度分布の表示
-
ホントは正規分布のハズなんだけど，
100 個ぐらいのデータでは上のようにぼこぼこ密度関数が憶測されてしまう．
こんなアヤしげな分布から
「直接乱数を発生させ」
たってインチキくさいだけなんだけど，
パズルとしてはおもしろい
(そして MCMC 計算結果の事後分布表示に多用してしまっている
density()
の性質に対する理解が深まるので)
ちょっと検討してみることにする．
-
そもそもこういうコトを検討する第一歩は
help(density)
(もしくは
?density)
しろってコトなんだけど，
そもそもこの質問者はマニュアルとかちゃんと読まないヒトで．
ともあれこの R の help を読むだけで density()
の挙動はいろいろわかる
……
と言いたいところだけど，
これはあるていどはカーネル密度推定なるものを勉強したことがあるヒトでないと
「このバンド幅って何?」
となるかもしれない．
-
より詳しい解説は例によって例のごとく
MASS4 本 5.6 節 Density Estimation
に書かれていたりする．
-
このぎょーむ日誌にまた数式画像ファイルとかはりつけて解説すんのは
面倒なので
(おお，
ゐきぺでぃあ
にはちゃんとした数式とわかりやすい図が示されているなぁ)，
「上のような憶測ぼこぼこ密度関数がどのように生成されているのか?」
を R コードで示すことで説明にかえたい．
ここで自作の
カーネル関数
K()
とカーネル密度推定値の vector fx
を (上の作図 R コードのつづきで)
K <- function(xj) dnorm(dxy$x, xj, dxy$bw) # xj は data.x のある値，ですな
fx <- apply(sapply(data.x, K), 1, mean) # それを全ての data.x の値について「ならす」
と定義してみる
(default では density() のカーネル関数は
正規分布なので dnorm() を使っている)．
さっきの図のつづきで，
このカーネル密度推定値の曲線を赤でを上書きしてみよう．
lines(dxy$x, fx, col = "#ff4000")
[density() 計算の再現]
もとの
density()
推定値の上にぴったり重なっていますよー，
とゆー図．
つまり
density()
の内部でやってることを
K()
関数と
fx 計算で完ペキに再現できてますよー，
ゆーコト．
-
さてさて，
density()
のしくみはかなり簡単だと確認できたので，
ここから乱数を発生させるような関数を定義するのは
簡単にできる．
たとえば rdensity()
と名づけた関数の実装例は
……
rdensity <- function(n, x)
{
dxy <- density(x)
K <- function(xj) dnorm(dxy$x, xj, dxy$bw)
rnorm(n, mean = sample(x, n, replace = TRUE), sd = dxy$bw)
}
ようするに，
データから
density()
計算して，
カーネル関数定義して，
-
データ内の一個をランダムに選ぶ (
sample())
-
そこを中心にして SD = バンド幅となる正規乱数を発生させる
という簡単なことだ．
これはくだらないと思われるかもしれないけれど，
density()
が生成する確率密度関数の定義と完全に首尾一貫している乱数発生方法になっている．
で，
くだらないついでにかくのごとく発生させた乱数と
さらにそこからもういちど憶測させた確率密度を図示してみる．
[rdensity() 乱数]
青の
|
が
rdensity() から得られた 100
個のサンプル．
さらに青線はその標本から再憶測させた
density()．
ウソにウソを重ねているというか，
どんどん「真の分布」から乖離していくような
……
-
いやー，
やってみたけどやっぱ
density()
からのサンプリングって意味不明でインチキくさいよ
……
えーい，
何やら邪悪な目的のために使うんだろうなぁ，
と懸念しつつ返信メイルかき．
rdensity()
の関数定義だけのせる．
説明は「明日公開のぎょーむ日誌みろ」
ですませておくか．
...
さて，ご質問の件ですが，
これは「density() が何をどうやっているか」がわかれば簡単にわかる問題で
す．答えだけ示すと
rdensity <- function(n, x)
{
dxy <- density(x)
K <- function(xj) dnorm(dxy$x, xj, dxy$bw)
rnorm(n, mean = sample(x, n, replace = TRUE), sd = dxy$bw)
}
これで rdensity(100, データの vector) などとすると御所望の乱数が 100 個
ばかり得られます．もうちょいくわしい説明は明日公開のぎょーむ日誌でも見
てください．この関数は density() の挙動に完全にあわせたモノになっている
のできわめて「正確な」乱数発生関数になっています …… ただし，私として
は density() なんぞで憶測されたアヤしげな確率密度関数を使って何かさらに
計算するってのがひどく危険なことのように思えます．とても邪悪な目的に使
われそうな気がしますね……
> コネタですが、某国立環境研PD氏にwinBUGSを布教しました。
> 今trap窓と激しい戦いを繰り広げているようです。
> 徐々に広がるべいぢあん・・・
貴君とヤツほど「ゐんばぐす Trap 窓地獄のどん底で子豚のごとくきーきーと
暴れまわる」ことが似つかわしい若手研究者もなかなかいないよなぁと思わず
にはいられません．
よいお年を．
さてススはらいは一件すんだ．
もひとつは昨日だかに書いた
「著者に何をいってもムダなので断りたい査読報告かき」
だが
……
1/8 提出の生態学会大会ポスター発表の要旨をひねくってみる．
2145 研究室発．
さすがに寒い．
買いもの．
2205 帰宅．
晩飯の準備．
晩飯．
[今日の運動]
[今日の食卓]
- 朝 (1020):
米麦 0.8 合．
ヒジキごはん．
ダイコン・ニンジン・ジャガイモの味噌汁．
- 昼 (1340):
米麦 0.6 合．
ダイコン・ニンジン・ジャガイモの味噌汁．
ゆでホウレンソウ．
ダイコンの酢のもの．
プチトマト．
- 晩 (2240):
米麦 0.8 合．
ダイコン・ニンジン・ジャガイモの味噌汁．
コマイ．
2007 年 12 月 31 日 (月)
-
1040 起床．
生活周期が破綻してるなぁ
……
朝飯．
コーヒー．
-
今年のススはらいをすべく，
例の断ってしまいたい査読報告処分を
……
と「なぜこの査読報告書きを断りたいか」
作文を書いてるうちにしんどくなって，
この作文をちょっと変えてそのまま査読報告にすればいいや，
と気づいた．
しかも前回却下時の (査読者ではなく) 編集者コメントがなかなか
おもしろい点をツいていることに気づいたので．
-
昼飯くって
1430 自宅発．
曇．
1445 研究室着．
-
査読報告かきに必要な統計学文献よみ．
-
大晦日においこまれてる院生の年越しこんさるを回避すべく
1800 研究室発．
1815 帰宅．
-
査読報告かきのつづき
……
なんだけど，
さすがにイヤになってきて進捗せず．
-
[今日の運動]
-
[今日の食卓]
- 朝 (1050):
シリアル．
- 昼 (1340):
シリアル．
リンゴ．
ヨーグルト．
- 晩 (2030):
北 13 条西 4 丁目
さか田
で買ってきた蕎麦
(500 円)．
