ぎょーむ日誌 2005-08-(11-20)

2005 年 08 月 11 日 (木)

• 0700 起床． ねむい． 昨晩は 2530 ごろ寝たんで． 朝飯． コーヒー． 0850 自宅発． 晴． 0905 研究室着．
• 昨日 の変態的ぎょーむ日誌を読んでアタマの中を整理すんのに時間とられる．
• これまた昨日， 低温研の加藤京子さんに母子里データ送っていただいたんで， 今日はそれで遊ぶことに． アタマがモデリングぐるいぎみだったので， 正気に回帰するにはよいかも．
• とりあえず， 舟尾さんの R tips の scatterplot3d・lattice 参考にしながら R で三次元図を描いてみる． まずはお手軽に scatterplot3d() (in scatterplot3d package)．

[scatterplot3d()]

scatterplot3d(
  d$X, d$Y, d$Z,
  pch = 20,
  color = topo.colors(max(log.rp))[log.rp],
  mar = rep(.5, 4)
)

おてがるなのはいいんだけど 「視点」の変更できない． Document にもそう明記されてる:

... but unfortunately it results also in a static plot, i.e. rotation is not possible.

• で， R で使える自由度の高い (視点変更など可能，など) 三次元作図関数群とゆーと lattice package なンだが …… これってかなり研究しないとうまい図が作れない． くそう Trellis graphics め．

[しくじり cloud()]

lattice package の cloud() 関数によるもの．

cloud(
  Z ~ X * Y,
  data = d,
  groups = log.rp,
  aspect = c(2, 3),
  points = Rows(
    trellis.par.get("superpose.symbol"),
    log.rp
  )
)

これでは色とかめちゃくちゃだ． なかなか難しい．

• 昼飯． 本日は札幌も暑い． といっても 30°C をちょっとこえただけなんだけど．
• 生態がっかいさーばー雑用． 久しぶりの (私にとっては) ちょー難問で 3 時間以上を費してしまった …… そしてまたかとーさんの助言にたすけられた． そしてますます TurboLinux Appliance Server という牢獄がイヤになってしまった．
• 今回の主敵は cgiwrap である． これは web server software (Apache) と連動して動くフィルターである． ネット経由で動かしたい CGI プログラムの owner がへっぽこだったりすると， 「お前のスクリプトなんか実行してやるもんか， ばーかばーか」 とばかにしやがる …… というわけでもないが， ともかくエラー吐いてとまる悪者である (ように見えるプログラムである)．
• しかも腐れ Turbo で構築した腐れ apache ではこの cgiwrap を必ず呼びだすよーに コンパイルしやがってるわけで …… つまり cgiwrap が切れないのである．
• いろいろ試行錯誤右往左往暗中模索七転八倒したあげくに， とりあえず当該ユーザーの権限を少し強めてやることで cgiwrap にばかにされない地位を与えた． これで CGI プログラムは動くようになった …… 以上の理解はどこか根本的にまちがった部分を含んでおり (以上のような対処でとりあえず CGI プログラムを動かせるよーになったのは 事実であるが)， もうちょいこの cgiwrap の動作を研究する必要がある．
• そもそもこの cgiwrap は私の報告にあるような人間の尊厳をそこなう邪悪な目的のためではなく， 本来は (計算機内の) 地位が低いへっぽこユーザーでも CGI プログラム動かせるように， という善意の目的で作られたハズなのだが …… あ， cgiwrap サイトに CGIWrap - User Instructions とかあるな． なんかすごく面倒なこと指示してるぞ．
• ばてた． とにかく自分で選べるときには， ぜったいに TurboLinux 系は使わん， と決めた．
• 小林さんからメイルいただく． 母子里ハナシではなく， アカマツ原稿 (光合成補償のほう) の修正に関して．
• しかし …… たーぼばてのせいか， こちらは進捗せづ． 1925 研究室発． 1950 帰宅． 晩飯．
• [今日の運動]
  • うんどう休養日
• [今日の食卓]
  • 朝 (0810): 米麦 0.7 合． キャベツ・ニラ・ホタテの炒めもの．
  • 昼 (1310): 研究室お茶部屋． 米麦 0.5 合． キャベツ・ニラ・ホタテの炒めもの．
  • 晩 (2230): 米麦 0.7 合． ネギ卵炒飯． モズク酢． 麻婆豆腐．

2005 年 08 月 12 日 (金)

• 0810 起床． 朝飯． コーヒー． 洗濯． 0925 自宅発． 晴． 0940 研究室着． けっこう暑い．
• とりあえずアカマツ原稿 (光合成補償のほう) の修正． メイル送信．
• 昨日， 低温研の加藤さんから再送していただいた母子里データ． R によるデータチェック & データ変換． ほーら， ちぇっくしてみたら「隠れた欠測点」 がぽろぽろと ……
• write.csv() はなんだか挙動がヘンであるように思うので write.table(..., sep = ",", ...) を使うことに．
• とにかく R の中でできる変換はすべてやってしまうことに． こんなかんぢで．

  > source("convCSV.R")
  ... (欠測点チェック) ...
  Warning: (NO LINE) NA at (3, 3.6, 9)
  Warning: (NO LINE) NA at (3, 4.2, 9)
  Warning: (NO LINE) NA at (3, 4.8, 9)
  # output to dataS.csv ...
  > (ds <- read.csv("dataS.csv"))
     Xmin Ymin Zmin Xmax Ymax Zmax IX IY IZ  L
  1   0.3  0.0  2.4  0.6  0.6    3  0  0  0 50
  2   0.3  0.6  2.4  0.6  1.2    3  0  1  0 54
  3   0.3  1.2  2.4  0.6  1.8    3  0  2  0 43
  4   0.3  1.8  2.4  0.6  2.4    3  0  3  0 47
  5   0.3  2.4  2.4  0.6  3.0    3  0  4  0 51
  6   0.3  3.0  2.4  0.6  3.6    3  0  5  0 61
  7   0.3  3.6  2.4  0.6  4.2    3  0  6  0 70
  8   0.3  4.2  2.4  0.6  4.8    3  0  7  0 62
  9   0.3  4.8  2.4  0.6  5.4    3  0  8  0 73
  10  0.3  5.4  2.4  0.6  6.0    3  0  9  0 41
  11  0.3  6.0  2.4  0.6  6.6    3  0 10  0 66
  12  0.6  0.0  2.4  0.9  0.6    3  1  0  0 46
  ...
  

• 昼飯．
• 1310 ごろより， またしても 「7F 『夜逃げ』退官教授」 たちが利用していた部屋のあとかたづけ作業 …… なんかそろそろ裁判所に告訴してやりたくなってきた． 本日は「図書室」の巨大本棚 2 セットの解体． 作業者 5 名． こういうばかばかしい作業はさっさと終了させるにかぎる …… と集中してとりくむ． 1520 作業終了．

[巨大本棚]

これは 2 セット目を解体中． サイズは 幅 0.6m × 奥ゆき 6m × 高さ 2.5m ぐらい． 私はボルト・ナットはずしに専従． 部屋が汚いのでマスク着用で作業． かなり汚れる．

• Gmail アカウント一個進呈． 残り 48 個．
• また母子里データ変換ぷろぐらみんぐ． 低温研の加藤さんにいろいろとお尋ねしつつ．
• データ構造がやはりたいへんで， しかも手ぬき明るさ計算のやりかたに依存する， と． このあたりうまく分離できんかな? 1830 研究室発． 1845 帰宅． 2000 自宅発北大構内走． 2055 帰宅． 体重 73.4kg． 晩飯．
• [今日の運動]
  • 北大構内走 55 分間．
• [今日の食卓]
  • 朝 (0810): 米麦 0.8 合． モズク酢． 麻婆豆腐．
  • 昼 (1240): 研究室お茶部屋． 米麦 0.5 合． 麻婆豆腐．
  • 晩 (2140): 米麦 0.8 合． ゴーヤ・ピーマン・エリンギ・ショウガ・イカの炒めもの．

2005 年 08 月 13 日 (土)

• 0730 起床． 朝飯． コーヒー． 洗濯． 怠業． 1020 自宅発． 晴． 1035 研究室着．
• いきなりプリンター問題． MacOS X 10.4 ``Tiger'' から CANON LBP-2810 にうまく印刷出力できません， とのこと． 調べてみると CANON の LIPS4 driver for MacOS X が Tiger に対応してないことが原因， と特定できた． ヴァージョンがちょっと変わっただけで printer driver の互換性うしなう Apple 社の OS もへっぽこだし， Tiger なんてずっと前にでた OS にまだ対応できない CANON も阿呆阿呆だ．
• またアカマツ (光合成補償) 原稿のちぇっく． 図とか作り直して ……
• とちゅうでわきみちにそれてしまった …… 先日 Gmail アカウント作る機会があったので， これを「オンラインスパムフィルター」 として使う設定と実験． Gmail でメイル転送やらせたいんだけど， スパム (いわゆる迷惑メイル) は自動選別して (Gmail には最初からこの機能がある) 転送しない， というふうにした． BananaBlog の説明など参考にしつつ ……
  1. (もしすでに転送設定してたら) いったん転送設定は解除
  2. 「設定」 でフィルタの設定に入る
  3. 「含めないキーワード」 に in:spam と書く
  4. 「次のステップボタン」 押す (警告でるけど無視無視)
  5. 「次のアドレスに転送」 で転送先を設定
  これでスパムが Gmail 「迷惑メール」内に残り， ふつーのメイルは転送される．
• とゆーことで， 私の「縦深メイル防御」システムにきわめて有効な一ステップが追加された． 現状はこんなかんぢで
  1. 北大あてにきたメイルはただちに Gmail に転送 --- 北大のメイルいんふらは信用できないので
  2. Gmail スパムフィルターでスパムを選別廃棄してから asahi-net に転送 --- このときコピーを Gmail に archive; あとで検索するときに便利だろうから
  3. いったん asahi-net に蓄積 --- asahi-net 使う理由は 「どこからでも送信」 (POP before SMTP) するのに便利だから
  4. 自分の端末から fetchmail で asahi-net のメイルとる (コピーは asahi-net に残さない) --- 蛇足ながら私が使ってる OS は Vine Linux である
  5. fetchmail で簡単なスパム防御 --- SMTP error handling
  6. fetchmail から procmail 経由で spamassassin に渡してスパム処理 --- これは Bayesian filter 型のもの
  7. 最終防衛線: 自作の mail viewer & handler とでもいうべき Perl スクリプトで Maildir の新着ファイルを「さっ」と表示して捨てる・残すの判断・処理 --- このへんがすごく「軽く」できるように作ってある; スパムでなくても即捨てできるメイルとかもけっこうあるんで
  8. メイラーで開く --- ホントに読んだり返信する必要のあるメイルだけ
  ともかく 2. の処理をいれることで， それ以降のスパム対策はほとんど不要になってしまった． いやー， すばらしい …… なにしろ私あてにくるメイルの九割超はすぱむなので．
• あ， よく見たら asahi-net も スパムブロック を無料で提供してるな …… まあ， Gmail のほうが検索なんかは強力だからそちらでフィルタリング & ためこみをやって， asahi-net は上のよーな使いかたに特化させるか．
• 後記， とゆーか …… 結局， asahi-net のスパムブロック設定することにした． というのも， こちらのアカウント経由でスパムが流れこんでくる， とわかったので．
• 北大がそういうすぱむ処理ちゃんとやれよな， と言いたくなるわけだけど (いちおー blacklist 方式のフィルタリングはやってるけど， これはかなりザル) …… いやはや， しかしながら， これまでの輝ける「実績」から憶測するに また何かとんでもない通信事故ひきおこすだけだろーから やらないほうがいいんだろな．
• 昼飯． 露崎さんがアラスカからもどる． あいかーらずカップ麺など食べておられる． 「起学」 (という名前の専攻科) の受験者 57 名ですか． 露崎さんいわく， 改組直後は「改組効果」で受験者が増えるものらしい． ちなみに， 私のいる (生物圏科学専攻の) 植物生態学コース は受験者 4 名だそーで． これは平年なみかな．
• アカマツ原稿修正のつづき． とりあえず修正文など小林さんに送る． 図 PDF ファイルのアップロード．
• で， アカマツ原稿全体のみなおし …… 作業の効率がすごーっく悪くて， ですね． 小林さんが書かれたところはまだしも匍匐前進的にススむとして， いつものことながらとゆーか自分が書いたとこにあたると， いかなる種類の逃避行動なのかよくわからんけど， いきなり脈絡なくむやみやたらと google 検索をはじめて止まらなくなったり wikipedia から流入する無用雑学をアタマの中に流しこんでみたり ゴキブリほいほいのぞくような気分で Gmail の 「迷惑メール」 フォルダ開いて 「おお，とれてるとれてる」 と感心したり ……
• は． 時刻はすでに 2000． なのに自分で書いた章のみなおしはまだ 1 ペイジもススんでないよ． 時間がふっ飛んでいるッ これはッ 何者かの スタンド攻撃 を受け ^{ゴゴゴゴ……} たわけないよな．
• ぢりぢりと進捗． こういう原稿には多々まちがいとかあっても不思議ではないんだけど， スタンド攻撃に気をとられてるせいか， ちっともみつからん …… よーやく発見， と思ったら私の説明不足に起因するものだったりして， ですね．
• とりあえず最後までチェック終了． 消耗． 報告はまた明日． 2225 研究室発． 2250 帰宅． 晩飯．
• [今日の運動]
  • いやはや休養日
• [今日の食卓]
  • 朝 (0830): 食パン． ゴーヤ・ピーマン・エリンギ・ショウガ・イカの炒めもの．
  • 昼 (1300): 研究室お茶部屋． 食パン． ゴーヤ・ピーマン・エリンギ・ショウガ・イカの炒めもの．
  • 晩 (2320): 食パン． ゴーヤ・ピーマン・エリンギ・ショウガ・イカの炒めもの． ニラ・ナメコ・豆腐の味噌汁．

2005 年 08 月 14 日 (日)

• 0740 起床． 朝飯． コーヒー． 怠業．
• 今日もいーぐあいに曇ってきた． 1100 自宅発北大構内走． 1225 帰宅． 体重 73.0kg． 昼飯．
• 1450 自宅発． 曇． 1505 研究室着．
• 院生たちとお茶部屋雑談． 明日からサロベツ実習 (これは起学の実習) でみなさんまたいなくなるのか．
• アカマツ補償光合成原稿， 昨日みつけたモデルまわり説明不足を修正する作文． ちょっとした作業なのにまたえらく時間かかる． メイルで送信． 時刻は 1625 か．
• なんだかちょっとアヤしげなところのある事務書類作成． といっても不正ではないんだけどね． ハンコをやたらと押す．
• 粕谷さん が統計学的手法の普及をいろいろと画策されているので， そのあたりのお返事かき． 粕谷さんは以前から， データ解析の夏 (もしくは春秋冬) の学校みたいなのが必要， というご意見だった (久保のご意見: 賛成だけど講師もうしつけられたら準備たいへんそう)．
• こういうのって学会大会とかのあとにやったほうがいいのかねえ …… 院生たちはけっこう旅費とかに苦慮してるわけで (上のあやしげ書類かきがそのあたりに関係あったりなかったり)． 大会みたいな「祭り」と無関係にやると， 集まれるヒトが少なくなりそう． 逆に講師のほうをあちこち転がすんなら 科研費 とかそのたぐいが必要だよね． そういう科研費のかてごりーってあったかな． 会議用のわくがあるかと思ったけどなかった． むりやりなら基盤 C とか?
• おお． Gmail アカウントまた一個進呈． 残り 47 個．
• そして今日も Gmail すぱむフィルターのおかげで たいへん快適．

[Gmail]

Web interface はこういうかんぢ． 他の webmail (asahi-net とか北大の) に比べて操作が快適という気がする (その代償として javascript つかいまくりなわけだが)． 検索もたぶんすごいんでしょう (まだメイルがそれほど蓄積してないんでよくわからない)．

• 母子里 MCMC 計算プログラミングは本日の進捗しなかった． 1940 研究室発． 2005 帰宅． 晩飯．
• [今日の運動]
  • 北大構内走 85 分間．
• [今日の食卓]
  • 朝 (0810): 食パン． タマネギ・ピーマン・卵の炒めもの． ニラ・ナメコ・豆腐の味噌汁．
  • 昼 (1255): うどん． タマネギ・ピーマン・卵の炒めもの． ニラ・ナメコ・豆腐の味噌汁．
  • 晩 (2150): 食パン． タマネギ・キャベツ・ニラ・ブナシメジ・ホタテの炒めもの．

2005 年 08 月 15 日 (月)

• 0800 起床． 朝飯． コーヒー． 0920 自宅発． 晴． 0935 研究室着． なかなか暑い．
• いきなり関係ない問題を検討しはじめる (じつは登校前からひたすら考えつづけてる)． ちょうど三年ほど前 に三次元空間内での 「円筒と直線の衝突」 とゆー計算問題を検討したわけだが …… これの計算量もうちょい減らせんかと考えてみた． とりあえずの別解としては
  1. 二直線間の最小距離 L と最近接点 X を求める (2×2 行列計算)
  2. 1. の最小距離 L が円筒直径以下なら 3. にすすむ
  3. 円筒両端中心点と点 X の距離が どちらも円筒軸長より短ければ「あたり」
  ただし， これはじつは円筒ではなく両端が半球になってる円筒みたいなカタチ， ということになってしまう (3. の計算において)． これを改良するには ……

[円筒と最近接点]

最近接点が円筒内にあるかどうかの判定には少なくとも二段階の判定が必要． 過程 3. の改良案として， 円筒端中心点 - 点 X - 円筒端中心点のなす角度 (の余弦) を余弦定理で計算し， これが円筒形状から決まる下限より大きい， とする方法もありえるな．

• アタマの準備運動がおわったことにして， 母子里問題にとりくむ．
• とりあえず Canopy クラスにパラメーターファイルとハコデータを読みこませて， 1300 個ばかりの Fbox を生成 (new) させるところまでできた． 昼飯．
• 昼飯後もプログラミングの続き． 三次元空間内に生成したハコたちを reference で結合していく． つなぎかたにも三通りあって， top_fbox (直上), effectors (光計算に影響するハコ), neighbors (葉が「移動」しうるハコ) ……
• 観測してる林分の一番上のハコの top_fbox はどうなるかとゆーと， これは dummy_top_box とゆーのを生成してくっつけている． このあたり欠測点のとりあつかいに注意が必要で ……
• あ， ここで dummy_top_box とやらでの光強度を 100% ではなく 99.9% に設定しとくと， 局所明るさの尤度評価のときにあれこれ都合よろしいな， という二項分布の性質を悪用した巧妙なごまかしワザを発見する．
• お茶部屋で一休み． 院生密度いよいよ低い． 絶滅危惧なのか? しかし試みに北大生協で「誘引エサ」買ってきて お茶部屋においてみると， 稀少野生院生たちが次々ととらっぷされていく． ということでお茶部屋雑談 …… 蛇足ながら， 小菅せんせーその他数名は夜行性なので昼間には捕獲できない．
• いんちき光計算を実現する世界の構築のつづき． 結合した fbox の重みの定義．

[方向と重み]

光量のもっとも単純なモデルを考えた場合， 面積あたりのエネルギーは sinφ に比例する (φ は「視線」の仰角)． かかる格子状配列の世界だと sinφ = sqrt(dz2 / (dx2 + dy2 + dz2)) と簡単化される (d* は視点座標と視線上点の座標の差)．

• にせデータを生成して， このへんの各種試験運転． 三次元空間ばぐがぼろぼろと見つかる．
• 修正作業してるあいだにいろいろと． 統計こんさるとか． 近ごろ卒業生たちのあいだで流行してる， ssh & りもーとぶらうざわざ， とか．
• 見つかったばぐは全部修正． とりあえず， 空間構造をあらわす 「配管配線工事」， 半分ぐらいは片づいたか． 2035 研究室発． 2050 帰宅． 2205 自宅発北大構内走． 2250 帰宅． 体重 73.4kg． 晩飯．
• また一個 Gmail 発送作業． どんどん Gmail 使いが増えてくね． ところで本日， 私の web browser から Gmail にうまくアクセスできんかった． ログインすると 「読み込み中」 とかで止まる． 一日中原因がよくわからなかったんだけど …… さっきふと思いついて browser のキャッシュを全部クリア したら， 問題なく Gmail にアクセスできるようになった． なかなかナゾ．
• 北大構内走中に上記「配管配線工事」の意外なめんどくささに気づく． てきとーにやっても問題なさそうに思ってたんだけど， じつはてきとーすぎると境界でヘンなことになったりするかも， という可能性ありそう． たとえば「はしっこ」で葉っぱが すかすかになるといった artifact だ．
• これを避けるためには葉っぱ 「移動」 なんかのときに 慎重な加重計算ともなう乱数発生が必要になるか． いやはや …… と書いたところで少しマシな対策を思いつく． 移動先の「加重」を数値ではなく reference の個数で表現する， という方式はどうかな? たとえば fbox1 の加重が 1 で fbox2 の加重が 2 なら 「行き先候補」配列を (fbox1, fbox2, fbox2) としておく， という方式． 境界めんどうも reference つけかえできりぬける， と．
• ちょっとそのあたり実装してみる． うまくいきそうだ．
• [今日の運動]
  • 北大構内走 45 分間．
• [今日の食卓]
  • 朝 (0850): 米麦 0.7 合． タマネギ・キャベツ・ニラ・ブナシメジ・ホタテの炒めもの．
  • 昼 (1240): 研究室お茶部屋． 米麦 0.5 合． タマネギ・キャベツ・ニラ・ブナシメジ・ホタテの炒めもの．
  • 晩 (2320): 米麦 0.7 合． ミニトマト． タマネギ・キャベツ・ニラ・ブナシメジ・ホタテの炒めもの． 麻婆豆腐．

2005 年 08 月 16 日 (火)

• 0700 起床． ねむい． 朝飯． コーヒー． 0900 自宅発． 晴． 0915 研究室着．
• 朝一番のかとーオフィスは暑い． 室温 32°C 超 (外気温は 23-24°C)． 理由は窓が東にあるから， そして計算機が多いからである． 私のまわりには二台 (web server の hosho と DNS server の lex) だけなんだけど …… かとー先生のまわりにはすごく古典的な Mac たちが 「成層」 して動いていて， ですね．

[堆積 Mac]

こういうかんぢで． しかもどのマシンも NetBSD といったハッカー御用達の OS なんかがインストールされていて， 稼働しているのである． この積もりかた， そして手前の SCSI (すかぢー) HDD なんかがいかにもはっかー， ってかんぢの風景ですね．

• まあ， まず前夜帰宅時にしめた窓あけて外気とりこんで， (さらに正午ちかづけば「日陰」になるんで) この部屋も 8F の中では涼しい一画になるんだけど．
• また北大の (学内→学外) 強制通過型くされ proxy server が腐れぎみだ．
• Gmail アカウント， また一個発送． Gmail アカウント作った人たちもしばらくすると， 「Gmail アカウント 50 個ばらまき権利」が発生するので， 世の中 Gmail 使いばかりになってよさそうなもんだけど， そうはなってないようである． ともかく使ってみたいヒトは連絡くだされ．

[Gmail account ネズミ講]

ところでさー， ひとつ不思議なんだけど …… 「紹介できるユーザーは残り 49 人です」 ってどゆこと? 私の計算では残り 45 のはずなんだが …… super spreader とみなされるとボーナスアカウント get? えーい， こうなったらどんどんばらまくか．

• 母子里 Gibbs 林冠計算のつづき． 計算空間のいんふらは整いつつあるので， とりあえずいんちきな初期状態を生成する過程にとりくんでみる．
• なンかヘンだな， と思ったら …… メートルとセンチメートルをとりちがえてるぢゃないか． はい， 今回の計算はぜんぶセンチメートルです， と観測データ変換過程にそういう計算追加する．
• Leaf クラスのコンストラクターを書いて， と．
• 1150 ごろ地震． A 棟 8F だと震度 3 ぐらいの揺れ? 腐れ proxy のせいで地震情報わからん …… あ， よーやくでた． 1146 宮城県沖地震， 仙台震度 6 弱， か． 東北大のヒトたちはたいへんだな．
• 昼飯．
• 葉群の初期配置生成ができてるようなんで， 図で確認するべく描画出力用関数を準備していく． 1445 できた．

[なんとなくな図]

ワクとかないと何がなんだかわからんな． そして葉群の生成に以外と時間かかることがわかった たぶん， Leaf おぶぢぇくと生成するときのメモリ確保が負担になってる． いや， 乱数による加重つきハコえらびに時間かかってるのかな．

• わくをつけてみる． これは座標系の問題あって多少めんどう． 作図やりなおすと， 今度は少しわかったよーな気分になれる． これは「ササ除去区」と呼ばれるダケカンバ林分の葉群 (の計算用の初期配置の一例)．

• 下は「ササ区」と呼ばれる林分での観測データを使って 生成した葉群初期配置の一例．

[葉群初期配置生成]

だいたい 600cm × 300cm × 700cm ぐらいのワクの中の 1300 箇所ほどで明るさが観測されている． 総葉数をてきとーにきめてやって (この例だと 10000 枚ぐらい)， 「明るさの上下落差」が大きいところに葉っぱが存在する確率は高い， というルール (重みつきランダム) で配置してやると，こうなる．

• これらはまだ「初期配置」， つまり出発点にすぎない． こういった配置を初期状態として Metropolis-Hastings algorithm (M-H 法) 適用する Markov Chain Monte Carlo (MCMC) 計算をススめていくことで， 観測データをよりよく説明できるような， さらに 「もっともらしい」 葉群配置たちを見つけだせるか (そのような葉群たちをどんどん生成していく Gibbs 過程の定常状態を発見できるか)， とゆーあたりがネラいだ． まだまだ先は長い．
• しかし三次元ばぐとりでアタマ消耗したんで休憩．
• ちょうどタイミングよく， とゆーか …… かとー家の片づけ手伝いの要請が． かとーさんが車でむかえにきてくれたので森君とそれに乗って， 北区某所のかとーハウスへ． 任務は 「ピアノの下のじゅうたんを交換したいんで， それを実現するべくピアノを (部分的に) もちあげる」 というもの．
• ピアノの先生である奥様の教室にもなる居間の三割以上を占居する グランドピアノなんだが …… これ， 重いね． 三本ある脚のうちひとつを持ち上げることはできるけれど， 三人がかりでも全体を宙にうかせることはできんかった． かかる制約状況のもと， ぱずるのごときじゅうたん交換作業を． とはいえ， 実際のところ 30 分ばかりで終了した． なかなか得がたい体験であった．
• 本職のピアノ運搬屋さんたちは， 専用のベルトと高度な技法もちいて二人でこのグランドピアノを 自由自在に移動させるらしい． すげー

[かとーハウス]

何か撮っても問題なさそうなモノ， とゆーことで …… かとーハウスの表札でも． なかなか立派な家でした． とうぜんのことながらかとーオフィスのように 散らかってるわけではない．

• 研究室までまた送っていただく． ばてたので， とゆー口実で院生たちとお茶部屋雑談．
• プログラムの細かい手なおし． あまり進捗せづ． 2110 研究室発． 2140 帰宅． 2205 自宅発北大構内走． 2305 帰宅． 体重 73.2kg． 晩飯．
• [今日の運動]
  • 北大構内走 60 分間．
• [今日の食卓]
  • 朝 (0800): 米麦 0.7 合． モズク酢．
  • 昼 (1255): 研究室お茶部屋． 米麦 0.5 合． 麻婆豆腐．
  • 晩 (2350): 米麦 0.7 合． ミニトマト． キャベツ・モヤシ・ネギ・ブナシメジ・サケあらの味噌汁．

2005 年 08 月 17 日 (水)

• 0820 起床． 朝飯． コーヒー． 洗濯． 0935 自宅発． 晴． 0950 研究室着．
• 昨日からてぬき光計算の具体的な手順について検討している． 一週間前 に示したように， 基本アイデアはこういうかんぢで．

• 方向と重みの関係は一昨日のこれ．

[方向と重み]

光量のもっとも単純なモデルを考えた場合， 面積あたりのエネルギーは sinφ に比例する (φ は「視線」の仰角)． かかる格子状配列の世界だと sinφ = sqrt(dz2 / (dx2 + dy2 + dz2)) と簡単化される (d* は視点座標と視線上点の座標の差)．

• 昨晩の北大構内走中， あるいは今朝おきぬけのぼーっとした状態で列挙した手順の断片あれこれ．
  • 明るさ減衰計算は全部 log(明るさ) でやる --- これによって正しく相乗平均的に計算できるだけでなく， Fbox ごとに計算をうまく分割できる
  • 各 Fbox が foliage vertical distribution (FVD) をもつ
  • FVD への加算・減算関数をつくる --- とりあえず枚数の増減だけがわかればよいはず --- いや， このあたり再検討が必要
  • FVD から三通りの対数減衰速度を計算しておく --- 加重ナシ (自分用)， 上加重 (他人用)， 下加重 (自分用)
  …… ありゃ． 意外とこれだけで MCMC 駆動に必要な てぬき明るさ計算は完結してるのかしらん?
• さて， こういう計算専用のクラスを作るかどうか， だが． うーん …… とりあえず LightProcessor クラス， とゆーのを作ってみるか．
• まずてはじめに， Fbox に全部おしつけてた作業の一部を LightProcessor に分担させればいいだけなんだが …… なぜかばぐとりにちょっと手間どる． 仕事のひきつぎを確認． 時刻はすでに 1125．
• 次第にタテわり官僚主義になりつつある母子里 Gibbs 林冠計算プログラムの構造改革のつづき． じつは光計算にまつわる仕事は， LightProcessor の「上司」である Fbox をトバしてもぎょーむがススむ場合が多い (つまり下っぱの LightProcessor どうしで情報交換すればよい)， とわかってきたんで， そのように「風とーしのよい」しくみに変更していく． 変更 & 試験運転のくりかえし．
• てなことやってるうちに， 上の図のような局所明るさ計算に入りこむ． LightProcessor で処理すべき計算はぜんぶ書けた． 時刻は 1355． 昼飯にするか．
• 北大生協に昼飯調達にでる． 研究室にもどって昼飯．
• またプログラミングのつづき． 昼飯前の作業がうまくできてるのを確認する試験運転． 次にこの計算世界でいちばんエラい Canopy から LightProcessor を順番によびだして明るさ計算やらせるところを作る．
• ケガれ言語 Perl のおぶぢぇくと指向なデストラクターって 使ったことなかったんだけど， 今回は使ったほうがよいかも， ということでデストラクター試験コードを作ってみる．

#!/usr/bin/perl -w
#-------------------------------------------------------------
package Leaf;
use strict;
sub new # コンストラクター
{
	my ($class, $id) = @_;
	my $self = { 'id' => $id };
	bless $self, $class;
	printf "Leaf %i is constructed.\n", $self->{'id'};
	return $self
}
sub DESTROY # こう書くとデストラクター関数として処理される
{
	my ($self) = @_;
	printf "Leaf %i is destructed.\n\n", $self->{'id'};
}
#-------------------------------------------------------------
package main;
use strict;
my %hash_leaf;
$hash_leaf{'1'} = Leaf->new(1); # hash に放りこんで……
delete $hash_leaf{'1'}; # 削除 
my @array_leaf;
push @array_leaf, Leaf->new(2); # array に放りこんで……
shift @array_leaf; # 削除
my $leaf = Leaf->new(3); # プログラム終了時にきえる $leaf
__END__

実行結果はこういうあたりまえなもので．

Leaf 1 is constructed.
Leaf 1 is destructed.
Leaf 2 is constructed.
Leaf 2 is destructed.
Leaf 3 is constructed.
Leaf 3 is destructed.

• これを今つくってるほうのプログラムの Leaf でも追加する． 生成するときにコンストラクターで LightProcessor に ->Add_leaf($self) して， 捨てるときに同じく LightProcessor に ->Delete_leaf($self) する， と．
• これで葉っぱの増減にともなう処理を簡単化できるんだけど …… いやはや， 私が上のようにしといたことをうっかり忘れて， ヘンなコードを書いてしまいそうでおそろしい．
• といった失敗しないように注意コメントを Delete_leaf 関数につけてるうちに， さらに強まった my マヌケさ (use strict に表現してみました) を発見してしまった． ここで Delete_leaf を呼び出せる状況においてはデストラクターは発動しない． なぜかと言えば， LightProcessor 内でその Leaf への参照が保持されているからだ (Delete_leaf はその参照を削除する関数)．
• ついでながら …… Perl の garbage collection は mark & sweep 方式．
• だめだ． DESTROY の使いどころない． いやはや． やはり葉っぱ捨て処理のときにまず
  1. $light_processor->Delete_leaf($leaf) してから
  2. %{$canopy->{'foliage'}} (蛇足ながらこういうのが Perl のケガれ書法だっ) に格納されている最後の reference を消す
  と． ややこしくて， 自分でもド忘れしそう …… なんで， ぎょーむ日誌に詳細を記してみる．
• とりあえず葉っぱ増減にともなう問題はあとまわしにして， 局所明るさ計算の試験運転とばぐとり．
• 局所明るさに関する試験運転， とりあえずできた． よい結果 & よくわからん結果でる． よい結果は計算時間が短かった， ということ． 葉数 10000 枚ぐらいだと最悪でも 1 秒ぐらいで計算がかたづく． これはこの先にひかえるめんどう MCMC 計算をススめるのに助かる．
• よくわからん結果はこれ．

• 昨日の方法でもって 「てきとーな初期状態」 (上の図の左) を作って， とりあえずその局所明るさ評価をやってみたんだが， それが右の図． いくら「てきとー初期状態」でもこれは明るすぎるよな． 葉数 10000 枚とゆーか LAI が 4 ってのはそんなに小さい数ではないはずなんだが．
• 計算まちがいか， あるいは手ぬき計算方法が手ぬきしすぎている可能性がある． ちょっと調べなおすか．
• えーと， 一個のハコ (60cm × 30cm × 60cm) に直径 10cm の球状の葉っぱ押しこんでいく場合 …… LAI = 4 を実現するには 91 個ほどつめこむ必要あり． で， 計算してみるとさすがにそれだけつめこめば明るさは ゼロちかくまで減衰する， とわかった． まず根本的なところはまちがってない， と確認．
• なーんか単純なミス， という気もするんだが ……
• しかし本日もケガれ言語プログラミングでアタマが消耗してんで撤退． 1920 研究室発． 1935 帰宅．
• しかしやはりヘンてこ計算結果は気になるんで気合いっぱつ バグ狩りの解析 …… 30 分で原因特定できた． 計算まちがいではなくプログラムのバグだった．
• いやはや， 上のほうで「風とーしのよい」しくみうんぬんといった仕事分割やった， と書いたけどそのときに入ったエラーだった． Fbox 間の接続がヘンになっていた． かくのごとく， 構造改革は難しい．
• で， これで明るさ計算やらせてみると …… おー， われながら「うまい初期状態」 生成させる方法を発明しちゃったよーで． この段階ですでにかなり「あってる」よーに見える．

• えーと上の図はどちらも横軸が ダケカンバ林分内の局所明るさの観測値で， 縦軸がモデルの「初期状態」葉群配置の仮想林分で 実現している局所明るさ (左は常数軸，右は両対数軸)． 明るさは 0-1000 であらわされている． 対数軸でみると観測値 50 (相対明るさ 5%) より下の観測点はほとんどないんだけど， モデル初期状態ではさらに暗い場所が林内に たくさん存在することになっている．
• まあ， たぶん「LAI = 4 ぐらい? だったら葉数は 10000 ぐらい?」 などなどといった前提のどこかがまちがっていたんだろう． そしてこの葉群「初期配置」の方法には何の正当化もない． なぜかとゆーと …… 私の直観だけで導出された計算方法であるからだ． 科学というよりオカルトですな． 実際のところまじめに検討すれば いろいろな欠陥を見つけられる．
• そこで …… こういった不適切な前提だとか， あるいは一見もっともらしいけどじつは根拠ふりーで 決めてしまっている葉の位置なんかを， 統計学・統計物理学の基本わざの連鎖わざでもって改善していくのが Markov Chain Monte Carlo (MCMC) 計算なのである． とゆーか， そうであってほしい， ということで今回の問題に取り組んでいるわけだ．
• たとえば， 葉数に関しては平均 10000 のポアソン分布を (ベイズ統計学ふうにいうなら) 事前分布として与えたとしよう． そのような状況での MCMC サンプリングの結果として生成される事後分布は たとえばもっと低い平均値をもつ確率分布になる， といった「改良」が得られるかもしれない …… そういう研究をやってるわけで．
• いやー， 「発注者」 小林さんが明日から札幌に来られるんだけど (その後， このモデリング対象である母子里の森林に調査へ)， その前に MCMC 計算以外の難所は なんとかすべて片づいた， か． めでたしめでたし．
• 一件落着したんで 2030 自宅発北大構内走． かなりよろよろ． 2125 帰宅． ばてばて． 体重 72.8kg． 晩飯．
• 本日の Gmail アカウント出荷は 1 個． まだまだありますよー
• [今日の運動]
  • 北大構内走 55 分間．
• [今日の食卓]
  • 朝 (0840): 米麦 0.7 合． キャベツ・モヤシ・ネギ・ブナシメジ・サケあらの味噌汁．
  • 昼 (1440): 研究室お茶部屋． 北大生協ツナサンド．
  • 晩 (2230): 米麦 0.8 合． メカブ． ニラ・エノキダケ・卵の炒めもの． キャベツ・モヤシ・ネギ・ブナシメジ・サケあらの味噌汁．

2005 年 08 月 18 日 (木)

• 0750 起床． 朝飯． コーヒー． 0920 自宅発． 曇． 0935 研究室着．
• 香川の小林さんが今日からこちらに出張とゆーか， 滞在研究というか私の追いこみに来られるので， かとーオフィスの空き机のまわりを片づける． いやはや， これが …… この部屋， A 棟 8F のゴミ捨て場として機能してたんで， あれこれ捨てるべきものが蓄積されてて， ですね． それらを地階まで輸送するのにてまどる． ひととーり終ってみれば時刻はすでに 1100．
• とあるかたから問い合わせ (とゆーかこちらからの押し売りとゆーか) で 「地図上の個体の位置にマークを起き， その下に個体番号をつける」 ということを実現する R スクリプトを書く．

tmap <- function(data, xlim = NA, ylim = NA, pch = "+", ...)
{
	if (length(xlim) < 2) # 作図領域 x の範囲チェック
		xlim = range(d$x)
	if (length(ylim) < 2) # 作図領域 y の範囲チェック
		ylim = range(d$y)
	plot(
		data$x,
		data$y,
		type = "p",
		xlim = xlim,
		ylim = ylim,
		pch = pch, # 点をあらわす記号の種類
		col = "black" # 点の色
	)
	grid() # グリッドを描く
	text(
		data$x,
		data$y,
		labels = data$No.,
		pos = 1, # 文字の位置: 1: 下, 2: 左, 3: 上, 4: 右
		...
	)
}

• 個体データが R の中では data という名前の data.frame おぶぢぇくとに入っているとしよう． 上の関数定義を R によみこませて， たとえば

  tmap(data, xlim = c(-10, 10), ylim = c(-20, 0), col = "blue")
  

  を実行すると下のような作図結果が得られる (tmap(data) だけなら全領域を作図する)．

• まー， 個体が密集してるとこでは文字が見づらくなってるけど， そのあたりはまた xlim = ..., ylim = ... で領域を指定して作図やりなおすしかないかな． プログラムをもっと洗練されたものにする， とう手法もあるけど．
• ともかくこのようなものを作ったのは …… なんでも某研究機関ではこういう図が必要なときは， 個体の位置だけをぐらふそふとで描かせて， あとはデータの表をニラみながら 全部手がきで 個体番号をつけてるそーで (個体数は数百もあるのに!) …… もはや怒りをとおりこして， その絶望的な悲惨さにコトバを失った． ということで上のような R コード書いて押し売りした次第である．
• 院生が 「海外の大学にメイル送ってもつっかえされてくるんですけど」 と報告してくるんで， 調べてみると ……
• 北大メイルサーヴァーがまたしても 「人類共通の敵」 と認定されていた …… これで三度目だ (一度目 と二度目)． ともかく， spamcop.net のブラックリストに登録されてしまったんで (理由はおそらく実際に北大から大量の迷惑メイルとやらが 世界中にまきちらされているから)， これが解除されないかぎりこちらからのメイル送信に いろいろ支障が生じるだろう．

[spamcop.net 認定 spam site]

iscan1.sys.hokudai.ac.jp (IP アドレス 133.87.1.96)． ここ で問い合わせることができる． これに掲載されると， spamcop.net 利用してメイルフィルタリングしてる機関に メイル送れなくなる (掲載→応答まで時間差がある)．
いつものことながら， まぬけな北大ネット管理者たちは この状況に気づいていない． しょうがないんで， メイル送って通報してやったんだけど， 返信もなければネット上でのお知らせもナシ． いやはや， また「なかったこと」にするつもりだな．

• そこに書いてあるところによると

  If there are no reports of ongoing objectionable email from this system it will be delisted automatically in approximately 5 hours.

  だそーで， 北大がまたすぱむ発信やらなければ 5 時間ぐらいしたら ブラックリストから自動的に削除されるらしい． うーむ， だいじょうぶなのか? 北大の管理者がさっさと解除申請しろよ．
• 昼飯． 昨年とかに比べると今年は大雪山も知床も登山者がちょい少ないそーで． 北海道登山はもはやハヤらない?
• Gmail のすぱむフィルター， なかなかいいデキだと思うんだけど完ペキではない． 感覚的には 20-30 にひとつぐらいはスパムを通してしまう， というかんぢだ (現時点では)． 私の縦深メイル防御体制の後衛に陣どる spamassassin はかなり対すぱむ経験を積んでいるんで， ここでたいてい迎撃できるわけだが．
• 院生の Dell デスクトップ機が起動しなくなったそーで． なーんか， 挙動からすると motherboard がコワれたんじゃないかしらん． で， フタをあけようとするんだけど， やりかたがわからず苦闘する． ひっくりかえしてみてよーやくわかった．

[Dell 機 open]

タワー型匡体左側を下に机におく． つぎに匡体上面と底面にあるボタンを押しながら， 匡体右側を持ち上げるよーにすると， このようにぱかっと開く． なんかカッコいいよね．

• で， HDD まわりがおかしいみたいなんで， いったん HDD ケイブルをぬく． 起動試験すると (今まで動かなかった) BIOS が動作しはじめた (HDD をつないでないんで boot しないけど)． それぢゃ， ということで HDD ケイブルをしっかりつけなおすと， こんどは問題なく起動した． どうやら HDD ケイブルの接触不良 (?) だったよーで．
• 他にもこまごまとした作業に …… ところで小林さんはいまだ現れづ．
• 1530 よーやく母子里林分 MCMC 計算にまた手をつける． 尤度計算に関するいんふら工事． R だったら何も考えなくてもできる計算あれこれの準備が ……
• 事前分布がポアソンってヤメたほうがいいかもな． variance が小さすぎるような気がしてきた．
• とりあえず MCMC 計算の基礎となる尤度計算の関数を作る．
• 北大ネット管理者からメイル． やはり「なかったこと」にするらしい． こいつらの情報隠匿はけっこう悪質だ …… 1730 ごろブラックリストからは削除された． 北大からのスパム発信が数時間とまっていたんで．
• 夕方， 小林さんがみえる． さっそく 8F 闇ネットに接続して， 何やらいそがしそうに働いておられる． その後， 母子里相談少々．
• 1850 研究室発． 小林さんと晩飯． 2005 帰宅．
• MCMC 計算の一ステップを書いていく． 葉っぱの確率論的な出し入れまわり． これはちょい難しいところだ． とゆーのも， 状態遷移を Metropolis-Hastings アルゴリズムの計算で reject された場合， すべてをもとの状態にもどす必要があるからだ． 詳細な計算を担当する LightProcessor に Make_backup & Restore_backup という関数をもたせて， これに対応してみる．
• すでに夜もおそくなってきたんでプログラミングは中止して， 2330 自宅発北大構内走． 走っているとアタマの中が整理される． 2420 帰宅． 体重 73.6kg．
• [今日の運動]
  • 北大構内走 50 分間．
• [今日の食卓]
  • 朝 (0840): 米麦 0.7 合． ニラ・エノキダケ・卵の炒めもの． キャベツ・モヤシ・ネギ・ブナシメジ・サケあらの味噌汁．
  • 昼 (1300): 研究室お茶部屋． 米麦 0.5 合． ニラ・エノキダケ・卵の炒めもの．
  • 晩 (1910): ヤマダモンゴル 北 8 条でジンギスカン (今日はラムの肩ロース) と飯．

2005 年 08 月 19 日 (金)

• 0750 起床． 朝飯． コーヒー． 0910 自宅発． 曇． 0925 研究室着．
• 北大メイルサーヴァー， いまだに「人類共通の敵」あつかい， つまり spamcop.net の ブラックリスト に掲載されたまま． 昨晩いったん解除されたはずなんだが， また北大からのすぱむ大量放出が観測されてしまったんだな．
• 解除されるのは早くても 14 時間後? 情報隠匿な北大ネット管理者たちもこの事件を 「なかったこと」 にするつもりだったんだけど， もはや隠しきれなくなったと観念したらしく， 「ブラックリスト (spamcop) に登録されております」 なるお知らせをこっそりと掲載している (学外からは見えないところに)．
• MCMC 計算にともなう確率論的な葉っぱだしいれ機構つくりのつづき．
• 1130 ごろ， ひととーりできる． MCMC 計算疾走 …… おお， ちゃんと葉っぱの数がそれっぽいところに遷移していく …… しかし， 葉っぱの出し入れにともなう Make_backup & Restore_backup まわりにまだ少しバグが残存しているな． ややこしいとこなんで詳細な見直しを始めてみる．
• 途中で昼飯． 小林さんと北大生協へ．
• ネヴァダ州ラスヴェガスから請求書くる． といってもばくちのそれではなく， アカマツ原稿の敵国語 editing のやつ． 34000 円 (5.5 時間)． 思ってたよりかなり安い．
• ついでにあやしげな書類も提出．
• またバグとりのつづき． ちなみに今やってる MCMC 計算は 葉っぱ一枚の出入りごとに Metropolis-Hastings 判定を適用するもっとも原始的な (同時にもっとも確実な) 方法である． なにしろてぬき光計算がかなり高速に処理できるようになってるんで．
• 低温研の加藤京子さんきたんで， 三人で夜まで母子里相談． 母子里 MCMC 計算に関しては
  • 葉っぱは球形ではなく円盤にしたほうがよい --- そのほうが他の葉数推定によくにた葉数になるだろう (ただししばらくは計算手間はぶくため球形を使う)
  • 円盤葉っぱにした場合， その葉方向モデルに工夫の余地があるかもしれない --- 理想的には MCMC 計算でそれっぽい角度分布が 「自動的に」推定される， という状況なんだが ……
  といったあたりか．
• 1950 研究室発． 加藤・小林・久保で晩飯にでる． 2110 帰宅．
• MCMC 計算のばぐ取りのつづき． undefined object エラーの原因はようやくわかった． カラのハコで葉っぱを減らすことはできない， というあたり前のものだった．
• しかしこの葉数がどんどん崩壊するのは …… たぶん単純なばぐだろうな． しかし原因がわからん． アタマがふらふらになってきたので寝る …… 後記: これは翌朝， 起床 10 分後に解決した． 寝てるあいだに問題が理解できたのだろう．
• [今日の運動]
  • うんどう休養日
• [今日の食卓]
  • 朝 (0840): 米麦 0.7 合． コマツナあえもの． 卵焼．
  • 昼 (1240): 北大生協食堂． 雑穀米 S． 豚汁 S． 韓国風春雨炒め． ナムル． 420 円． 644 キロカロリー．
  • 晩 (2010):

[ばら寿司]

北 13 西 3 のサカナ系定食が良いごはんや デンスケ のばら寿司 (金曜日だけ，数量限定)． マグロ・サバ・カニ・イカ・イクラなどなど． うまい． 750 円． 安い．

2005 年 08 月 20 日 (土)

• 0700 起床． 雨． 起きていきなり母子里林冠 MCMC 計算のプログラムながめてたら， なぞの葉数減衰ばぐがわかった． 葉っぱの「戸籍簿」管理がしっかりしてなかった， だけだった． 朝飯． コーヒー．
• で， よーやくにして MCMC 計算がまともに疾走するようになって， いろいろなことがわかった． とりあえず計算結果の例． 葉群配置は「てきとーに作った」初期状態 (下の図左) と MCMC sample の例 (右) ではそれほど変わっているようには見えない． しかしながら ……

• とゆーか当然のことながら， 三次元区画内に格子状に配列された 1320 個の光センサーの 計測値は初期状態のばらつき多い状態から MCMC 計算によって 大幅に改善されている． 現場での観測値をかなりうまく再現したものになっている．

• Metropolis-Hastings (M-H) アルゴリズムを使った計算進捗状況はこんなかんぢで． たとえば， 初期状態に対する最初の改良 step1 においては 葉っぱ消し去り (Decrement) を試み， それによって対数尤度 (に比例する量) が -31537 から -31521 に改善され ACCEPT されて新しい状態として採択されている．

  # Decrement: (-31537 vs -31521) ACCEPT     1: ll = -31521.7, nl =  7562
  # Increment: (-31521 vs -31524) REJECT     2: ll = -31521.7, nl =  7562
  # Increment: (-31521 vs -31595) REJECT     3: ll = -31521.7, nl =  7562
  ……(中略)……
  # Increment: (-7706 vs -7741)   REJECT 99994: ll = -7706.7,  nl =  7092
  # Increment: (-7706 vs -7706)   ACCEPT 99995: ll = -7706.6,  nl =  7093
  # Decrement: (-7706 vs -7979)   REJECT 99996: ll = -7706.6,  nl =  7093
  # Increment: (-7706 vs -7707)   ACCEPT 99997: ll = -7707.2,  nl =  7094
  # Increment: (-7707 vs -7709)   REJECT 99998: ll = -7707.2,  nl =  7094
  # Increment: (-7707 vs -7715)   REJECT 99999: ll = -7707.2,  nl =  7094
  

  いっぽうステップ数が 10 万にちかいこのあたりでは， 対数尤度の変化はほとんどなくなり， Gibbs 過程の定常状態といってもよい状態になっている (このあたりでは尤度は単調に改善されず， M-H 法ではある確率で「悪化」することも許容されている)． このような状態の林冠たちを sampling していく (そして林冠の光合成生産のみつもりなどに使う)． 定常状態における葉数は 7090 枚あたり， となった (初期状態では 7563 枚だった)．
• わかったことあれこれ．
  • 葉っぱ増減のやりかた: 「こうしたら効率いいだろう」 と考えてたのはそれほど良くない， とわかった --- ランダムにハコを選んででたらめに葉っぱ出し入れやって， M-H 法の ACCEPT/REJECT 裁定にまかせる， といったしんぷるな手順が良いようだ
  • 計算にすごく時間かかる --- このように葉っぱ 7000-8000 枚ぐらいだと， 定常に到達するのに 10 万ステップぐらい必要?
  • 初期状態やはり重要 --- 言うところの「てきとーに決めた」初期状態は (他のさらにでたらめな初期値に比べると) 格段に尤度が高く， これのおかげで定常に到達するまでの時間は節約できてるよーな 気がする．
  • 葉数などの (明示的な) 事前分布は必要か? --- この問題ではポアソン分布などで 葉数の事前分布を仮定しないほうが良いのかもしれない …… 結果がこれにかなり左右されるんで (ということで， いま使ってるのは無情報事前分布)
  • 収束判定は難しい --- Full model や Null model の尤度を計算して， deviance みたいなのも収束判定の参考にすべきかもしれない
  いやー， 面白いけどたいへんですなぁ．
• とりあえず一件落着したんで 1030 自宅発北大構内走． 曇天ときどきすこし雨． 1150 帰宅． 体重 72.2kg． 昼飯．
• 霜月さん ところで少しらくがき． R で組みあわせ列挙するなら， たとえば

  library(gtools)
  combinations(5, 3) # 5 個の中から 3 個とる組み合わせすべて
  

  てなかんぢで (CRAN の gtools が必要)， ケガれ言語 Perl のケガれ書法で同じこと書くなら

  use Math::Combinatorics;
  print join("\n", map { join(" ", @$_) } combine(3, (1 .. 5)));
  

  としてみたり (CPAN の Math::Combinatorics 必要)．
• 1300 自宅発． 雨． 1315 研究室着． 1325 研究室発． ますますどしゃぶり状態． 北大クラーク会館地下の 理髪店 で散髪． 1430 終了． 2650 円． 雨あがっている． 研究室もどる．
• 研究室の共用マシン Dell 機 (Pentium4 2.4GHz) で上の MCMC 計算やらせてみた． 私の ThinkPad X31 (PentiumM 1.6GHz) だと 55 分かかる 10 万ステップ計算が 22 分で終了した． 速い．
• 計算の顛末などをかくのごとくうだうだとぎょーむ日誌に書いてみたり．
• お茶部屋で井田君と モデル選択解説論文 (Johnson and Omland, 2004) の話とか …… だー， また CANON のへっぽこ印刷機が トナーぎれか． ひとつの色がなくなっただけで印刷不能になるってのもなぁ．
• 母子里林冠 MCMC 計算について理解がススんだんで， それにあわせて計算プログラム (in Perl) りすとらしてみる． よけーな変数や関数を削除していく．
• 母子里とは無関係の， 上の Combinatorics (combination 列挙) 問題にちょっとはまりこんでしまう．
• 2030 研究室発． 小林さんと晩飯． 2140 帰宅．
• Perl 自作関数による combination 列挙問題の回答例つくって， 霜月さん ところに投稿． えー， 問題は N 個の中から k 個とるすべての組み合わせを列挙せよ， というもの． たとえば N = 5 で k = 3 だと

  # 実行結果
  1, 2, 3
  1, 2, 4
  1, 2, 5
  1, 3, 4
  1, 3, 5
  1, 4, 5
  2, 3, 4
  2, 3, 5
  2, 4, 5
  3, 4, 5
  

  と示せば良い．
• 解決策で一番簡単なのは， 上で述べたように R だろうが Perl だろうが 利他的なヒトたちのつくった library を利用すること． しかしながら， それがイヤだという場合はどうすればよいか． Perl だとこういう combo なる関数を定義したらいいんじゃないの， とゆーことで ……

  #!/usr/bin/perl -w
  use strict;
  sub combo 
  {
  	my ($k, $array_orig) = @_;
  	my @array = @$array_orig; # copy
  	my @set;
  	if ($k == 1) {
  		@set = map { [$_] } @array;
  	}
  	else {
  		for (0 .. (scalar(@array) - $k)) {
  			my $head = shift @array;
  			push @set, map { [$head, @$_] }
  				@{&combo($k - 1, \@array)};
  		}
  	}
  	return \@set;
  }
  # 実行例 1
  print join("\n", map { join(", ", @$_) } @{&combo(3, [1 .. 5])}), "\n\n";
  # 実行例 2
  print join("\n", map { join(", ", @$_) } @{&combo(3, ["a", "b", "c", "d"])}), "\n\n";
  # 実行例 3
  my $i = 1;
  for my $s (@{&combo(4, [1 .. 8])}) {
  	print $i++, "\t", join(", ", @$s), "\n";
  }
  

• 考えかたの基本は再帰で (&combo(3, [1 .. 4]) の場合)

  	-+-1-+-2-+-3
  	 |   |   +-4 
  	 |   +-3---4
  	 +-2---3---4
  

  とゆー「木」を作ればいいのである．
• ぢつは似たような考えかたで解決する問題について， すでに 去年 1/4 とか 去年 8/17 あたりですでに考えていたもんで．
• [今日の運動]
  • 北大構内走 80 分間．
• [今日の食卓]
  • 朝 (0810): 米麦 0.7 合． コマツナあえもの． 卵焼．
  • 昼 (1240): うどん． コマツナあえもの．
  • 晩 (2100):

[冷麺]

北 12 西 3 の 喜和苑 の冷麺定食． 950 円． 麺は冷麺のそれではなくふつーなんで， 冷やし中華とゆーか． 味つけは「硬派」なかんぢでうまい． 自家製キムチいいんだけど， ちょい量が多い (と塩分とりすぎを気にしてみる)．