ぎょーむ日誌 2002-12-(11-20)

2002 年 12 月 11 日 (水)

• 0800 起床． 朝飯． コーヒー． 0850 自宅発． 雪． しかし路上の積雪はそれほど多くない． 気温はすごく低い． 先週末からずっと真冬日が続いている． 0900 研究室着．
• ここしばらく mlfitting に関するご指摘をいろいろと受けたので， 定番数値計算本 Numerical Recipes などを読んで次なる改良について検討する ……
• ということで午前中は Polytope 法の復習をして， 今までナゾだった修正 Powell 法 を新しく理解できたところで終わった． Numerical Recipes の example ってなんとも難解なる 書きかたなもんで ……
• 路上にうっすらと積雪しさらにぱらぱらと降雪つづく 北大構内に走りに出る． 1200 研究室発． ゆっくり走れば自分がコケる心配ない程度の路面状態なんだけど， 自動車がよろよろとこちらに走ってくるのは何ともおそろしい． 1235 もどる． 昼飯．
• 1300 からのはずが 1315 から 「この本を書いたときに Bazzaz 先生は どれぐらいくたびれきっていたのか 憶測する会」． 今日の担当はかとー先生． 敵国語はっきんぐも得意とする かとーさんなので文の解釈には問題ないのだけど …… Chapter 8. ``Competition and the evolution of response breadths and niches'' はこれまでの中で少なくとも (図と本文の不一致など) 単純ミスはもっともヒドく， 内容も支離滅裂ぎみと言ってよいだろう． Bazzaz 先生はともかくかとー先生はえらく消耗してました． まぁ， この本は著者の強引なるこじつけはてきとーに読み飛ばして， 北米放棄農地にわさわさと生えてくる 草本たちの不思議な挙動を楽しめばよいのである． 1600 終了．
• 昨日着手したゑくせるファイルなる virus を 自動解毒するプログラミングを継続する． 今日わかったのは元ファイルに日本語文字が混入している場合，

	use Spreadsheet::ParseExcel;
	use Spreadsheet::ParseExcel::FmtJapan;
	my $formatj = Spreadsheet::ParseExcel::FmtJapan->new(Code => 'euc'); 
	my $file = shift @ARGV;
	my $oBook = Spreadsheet::ParseExcel::Workbook->Parse($file, $formatj);

というふうに Spreadsheet::ParseExcel::FmtJapan モジュールを使う必要がある．
• さらに悪辣なる機種依存文字のたぐいで汚染されている場合は ...::FmtJapan2 を使わねばならんのだけど， なんだがこっちはエラー出て動かない． えーい， ...::FmtJapan では解毒作用が低くて「シュート」が「シュ〓ト」になる． 入力に使った M$IME とやらがあいかわらずへぼくて「ー」に おきてやぶりな文字コードも吐くんだよな． 過去の失敗反省せぬ阿呆企業め． まぁ， このあと別の Perl モジュールに読ませて データベイス化するわけだし 「シュ〓ト」 のままでいいや．
• 以上の日本語フォーマットをのぞけば， あとは CPAN サイトだので示されている Spreadsheet::ParseExcel の用例のように動かせる．
• わーくぶっくとやらをばらばらにしてテキストファイルとして 一枚一枚の内容を個別に取り出すコードは書けたんだけど …… やはりゑくせる使いはファイル作成も力まかせなので， 読みだしプログラム作るほうとしては苦労させられるわけだ． とはいえ， 以前に解毒せねばならなかった激烈に凶々しい物件あれこれに比べれば よほどマシだが．
• 不要な部分をばさばさと切り捨てるしくみを作っていくと， 次第にデータ格納フォーマット不統一が見つかりはじめる． ゑくせるファイルではありがちなことだ． このあたりはデータ持ち主の M1 大澤君に 直してもらわんことにはどうしようもない． ということで今日はここまで．
• 帰るついでに浦口さんカエデ計算の進展をいろいろと教えていただく． おお， 想定条件への応答が樹種によってかなり異なってみえる． 方法論的な試行錯誤はこれでだいたい片づいたんではなかろーか． ただし成長モデル推定などに関してはまだ改善の余地あるけど． 個体群という概念をあからさまには待ち出さぬ 確率論的な樹木の生涯， というのはどういう提示のやりかたがウケるだろうか …… 2045 研究室発． 北 12 生協で買いものして 2100 帰宅． 体重 73.2kg． 晩飯．
• ゑくせるファイル自動解体作業ばて．
• やっぱり昨日のあやしげなる数式は間違いだったそうで ……

	xy <- scan("sample.txt", list(x = 0, y = 0))
	require(nls)
	Pmax = 10.2 # きめうち
	psynthesis <- nls(
		y ~ (f * x + Pmax - sqrt((f * x + Pmax)^2 - 4 * (f * x * q * Pmax)))
			/ (2 * q), # さらに葉の呼吸も引くのかな?
		data = xy,
		start = list(f = 0.1, q = 0.5) # とりあえず
	)
	summary(psynthesis)

• [今日の運動]
  • 北大構内走 1200-1235． ストレッチング．
• [今日の食卓]
  • 朝 (0830): 米麦 0.8 合． ナバナ・タマネギ・マイタケ・卵の炒めもの．
  • 昼 (1245): 弁当． 研究室お茶部屋． 米麦 0.8 合． ナバナ・タマネギ・マイタケ・卵の炒めもの．
  • 晩 (2140): スパゲッティー． タマネギ・ニンニクの炒めもの．

2002 年 12 月 12 日 (木)

• 0710 起床． 朝飯． コーヒー． 0850 自宅発． 曇． 0900 研究室着．
• 昨日うまくいかなかった 悪辣なる機種依存文字汚染対策である Spreadsheet::ParseExcel::FmtJapan2 の問題点わかった． (私の場合だと) /usr/lib/perl5/site_perl/5.6.1/Spreadsheet/ParseExcel/ に FmtJapan2.pm がインストールされてるんだけど， これの 21 行目から

	my $oMap = Unicode::Map->new('CP932Excel');
	die "NO MAP FILE CP932Excel!!" 
		unless(-r Unicode::Map->mapping("CP932Excel"));

というコードがある． これら 3 行から ``Excel'' なる邪教のコトバを削除すると 正常に動作する …… 理由はやはり CPAN モジュールのひとつである Unicode::Map でインストールされる 解毒 map の名前が CP932Excel.map から CP932.map に変更されたため．
• 今日は夜までかかって， この「観測データをゑくせるファイル内から救出する」 プログラムのしあげをした．
  • 毎週の調査ごとにゑくせるファイルがひとつ存在する．
  • ひとつのファイル内には複数のシートが存在する．
  • ひとつのシート内には複数の表が存在する．
  • ひとつの表には最新のものだけでなく 前回・前々回のデータが書いてあったりなかったりする．
  というような二重三重のカベにはばまれた獄舎のごとき劣悪なる状況 (生態学分野における多くのゑくせりあんは かかる悪質で変態的な拉致監禁を好む病癖をもつ)， そこから必要なデータを助け出し (日付判定などを行っている)， ひとつひとつの表ごとに個別のテキストファイルとして保存する Perl スクリプトができた， ということだ．
• 20 余週間ぶんの呪われ.xls ファイル 20 余個をまとめて このスクリプトに放りこむと …… 3 分間ほどかけて 200 個以上のテキストファイルが生成された． ざっと見たところ 変換はおおむね正しいように見える． しかしデータベイスに格納する段階で不具合が発見されることは 間違いないから， そのときはまた元データにもどって直さないといけないだろう．
• しかしばてたので本日はここまで． 1830 研究室発． 雪がどかどか降ってる． また東京都から納税しろ命令がきたので 札幌市役所ちかくの ATM まで現金をおろしに行く． ついでにひさしぶりに本屋などふらふらとしてから 2000 帰宅． 体重 73.0kg．
• どうもばてたなぁ．
• [今日の素読]
  • David Salsburg. 2001. Preface, in ``The lady tasting tea''. Owl Books.
• [今日の運動]
  • 雪の北大構内走 1210-1245． ストレッチング．
• [今日の食卓]
  • 朝 (0820): 米麦 0.8 合． コマツナ・タマネギ・ネギ・卵の炒めもの．
  • 昼 (1255): 弁当を自宅に忘れてきたので 研究室お茶部屋でパンをかじる．
  • 晩 (2020): 昼に食わなかった弁当． 米麦 0.8 合． コマツナ・タマネギ・ネギ・卵の炒めもの．

2002 年 12 月 13 日 (金)

• 0820 起床． 朝飯． コーヒー． 0900 自宅発． 晴． 積雪 5cm ぐらいか? 雪面からの反射がまぶしい． 0910 研究室着．
• 新しく apt-get install した Canna 3.6p1 の設定． この Unix 系その他 OS 用の free かな漢字変換サーヴァー Canna では個人別頻度学習ファイルの使いまわしはできないんで ……
  1. for dic in `lsdic`; do rmdic $dic; done で古い頻度学習ファイルを捨てる
  2. for dic in `lsdic -i`; do mkdic -fq $dic; done で新しい頻度学習ファイルを作る
  3. cat usr1.ctd | mkdic - user 退避しておいたユーザー辞書を使って 新しいユーザー辞書の初期化
  4. /usr/doc/Canna*/README.* や /var/lib/canna/dic/canna/dics.dir など見ながら $HOME/.canna の辞書指定を変更 (iroha やめて gcanna, gcannaf にするとか)
  5. kinput2 を止めて再起動 (KTerm とかも再起動しなくては …… Rxvt は再起動しなくてよい?)
• 昨日はばててたので書けなかった Spreadsheet::ParseExcel 報告のぎょーむ日誌書き．
• D 論予備審を無事に終了しあとは旦那のいる欧州に行くばかりのはずの ユニさんが図版印刷出力もってきて「ちょっと教えてください」 とのこと． なんでもフォントじゃぎー嫌いの甲山さんが じゃぎーな図を全部さしかえろとか言ったらしくて …… で， その作図環境というのが (まず Mac 内で) ゑくせる→なにかグラフ作成ソフトウェア→ Canvas5 とかいう描画ソフトウェア→ (ファイルサーヴァーの Linux 機経由して) →わーど2000 for Windows． これまた生態学者にありがちな 複雑怪奇なる全手動な作図・組版工程である．
• Canvas5 から出力してるところでじゃぎーになってる． なぜか Mac や Windows のソフトウェアでは EPS ファイル取り込みというような ふつーの機能が欠落しているためだ． わーどの呪われコンヴァーターの挙動を眺めてるうちに， Mac 独自な PICT 利用に可能性を見いだす． Canvas5 に PICT 出力させてわーどに取り組むと …… かんぢんなところで手抜きしやがるんだよな． じゃぎーは解呪されたんだけど， 上つき文字列 (superscript) が二重描画，か． しょうがないので 「superscript 部分だけ消去して， その部分を小フォントの別文字列にしてください」 ということで． 作図工程を自動化してないとこういうとき面倒でしょう．
• 自由集会 の要旨を送る …… え? もっと短く? ということで 3 割ほど文言を削って再提出．
• 1210 雪の北大構内走発． 晴れ． いったん北大構内から出て銀行に行き， 東京シンタロウなる地頭への上納金 52000 円を支払う． もいっぺん支払えさもなくば軍団を派遣するぞと脅されてるので， 年に合計 20 万円以上を巻き上げられる計算になるわけか． おそるべし武蔵国の年貢． 北大構内をひとまわりして 1250 研究室帰着． 昼飯．
• どうも午後から仕事ススまぬ． 昨日， 救助ルートを確立したデータたちを 今度は Perl オブジェクトなデータベイスにしないといけないんだが …… prunusdb.pm， すでにホネぐみはできてるんだけどねぇ．
• ということで， 浦口さんの仕事をじゃまする …… カエデ成長モデルに関して， 何点かよくわかってなかったあたりの質問して教えていただく ……
• なンと驚愕すべきことに， 浦口さんは自分でベータ乱数を生成する Perl モジュールを ネット上で公開されてる C のコードをもとに書いていた ……! あとから伏見正則「乱数」本 (UP 応用数学選書) を見ると， Kinderman-Monahan-Ramage の算法， というやつのようで． いやー， すごい． Perl 作者の Larray Wall も仰天するだろう．
• この成長モデルに関して疑似尤度 (quasi-likelihood) 法をを用いて， 不都合な部分を回避してみようという策略の共謀． 満たすべき条件を確認してから
  • 疑似尤度法で 「分散は (たとえば) 平均値に比例」 という仮定のもとで平均に関連する パラメーターたちを推定する (確率分布は特定せずに)
  • 成長サブモデルを説明するときに 指数関数族から都合の良い確率分布を ひとつ特定して使うことにする
  …… という方針でおしてみることにする．
• Quasi-likelihood についてひととーり手元の本とか ネット上の説明をチェックしてみる． とりあえず， 対数線形モデルにそって非等分散正規乱数を生成してから， それを標本集団として R の glm で family = quasi(link = "log", variance = "mu") というふうに指定して推定計算させてみる．

	dataset <- scan("sample.txt", list(x = 0, y = 0))
	dataset <- list(
		x = dataset$x,
		y = dataset$y,
		logx = log(dataset$x),
		sqlogx = log(dataset$x)^2
	)

	quasimodel <- glm(
		y ~ 1 + logx + sqlogx,
		family = quasi(link = "log", variance = "mu"),
		control = glm.control(maxit = 20, trace = TRUE),
		data = dataset,
		start = c(-2.0, 1.0, 0.2)
	)

	summary(quasimodel)

• ところがいかんせん， この推定問題が難しすぎるためか R の根性が欠落してるためか， glm はこの状況では推定値を収束させることができないようだ．
• いろいろぢたばたしてみたんだがうまくいかない． 先日，森林総研関西の宮本さんが熱帯林現存量推定で R の glm を使ってうまくいったのは …… なるほど， 実験してみると平均値に関するモデルが 単調な増加関数だったりすると， わりとうまく推定してくれるみたいだなぁ． ところがあのわがままなるカエデの成長はふにゃふにゃしてるんだよねえ．
• 計算ばてしたので現況報告してから 2050 研究室発． ところが A 棟の階段を降りていくうちに quasi-score function がアタマの中で 「ほどけて」くる …… で 1F に到達したときには汎用最尤推定な mlfitting でこの問題が解決できるんでは， という見とおしが多少よくなった． 2100 帰宅． 体重 73.6kg． うーむ減らん．
• 晩飯の炊飯しておいて北 12 生協で買いもの． ついでに電話代 1900 円・電気代 1179 円の支払い． 晩飯．
• 晩飯後にでれでれしてから， 疑似尤度の正体をみきわめる数値計算など．

[Quasi-score 関数] てきとーに書いたもの． 分散∝平均モデル．

• よーし， これなら mlfitting の plugin_function.cc を書けば自分で計算できそうだ， ということで取り組んでみる．
• 2420 計算結果得られた． なかなか悪くない．

[Quasi-likelihood 推定] 非等分散正規分布 (分散 = 0.25 × 平均) な乱数集団への mlfitting による非線形モデルのあてはめ (quasi-score の最大化)． 青線が推定結果で緑線が真の平均．

• しかし分散の係数をどうやって計算するのが妥当なのか， よくわからない．
• [今日の運動]
  • 北大構内走 1210-1250． ストレッチング．
• [今日の食卓]
  • 朝 (0840): 米麦 0.8 合． ナバナ・タマネギ・マイタケ・豆腐の炒めもの．
  • 昼 (1300): 弁当． 研究室お茶部屋． 米麦 0.8 合． ナバナ・タマネギ・マイタケ・豆腐の炒めもの．
  • 晩 (2150): 米麦 0.8 合． ハクサイ・ニンジン・タマネギ・豆腐の味噌汁．

2002 年 12 月 14 日 (土)

• 0800 起床． 朝飯． コーヒー．
• また quasi-likelihood 推定法の計算の確認のつづきなど …… 1040 自宅発． 曇． 1050 研究室着． 今日は北大・低温研で生態学会北海道支部の 地区会 あるんでそちらに行く． 1100 研究室発． 1110 低温研着． 一般講演がおわり院生たちの発表が始まるところだった．
• 液晶プロジェクターのごたごたあれこれ． いやはや．
• 1700 すぎまで北海道内にいる院生たち発表 (ほとんど M2 のもの ……つまり修論発表会状態) つづく． 函館の水産学部でなかなか面白そうな研究やってる． カニの精子制限とオスの繁殖戦略とか 貝殻形状 = 遺伝 + 環境とか．
• 私のいる研究室の M2 諸氏の発表は …… いまいちだなぁ． それなりに準備はしたのはわかるけど．
• 彼らにかぎらずゆーい差決戦主義者が多い …… ゆーい差決戦主義とは， 統計学的ゆーい差→生態学的有意差なり， と解釈するカルト宗教． いわゆる「P 値」が「1 - 自分の言ってることが正しい確率」 あらわすと狂信してるので， この P なる数値が小数点以下何ケタであるか をひどく重視する． P < 0.000001 とかいう計算結果が得られようものなら かなり嬉しそうにそれを示す．
• しかし「ゆーい差ありました」 ぐらいしか言えない研究も多い． つまり観察した現象を記載してるだけ， というか． 記載といえば， 群集構造の「分類」で TWINSPAN (two way indicator species analysis) とやらがやたらと普及してるな …… それってまともな統計モデルにもとづいてるのか? ちゃんと計算ソフトウェアの動作とか チェックしてるんだろうか．
• 1700 すぎに発表が終わったので低温研から撤収． かとーオフィスにもどって quasi-likelihood 推定における ばらつきパラメーターの推定法の検討．
• あれこれ考えたり試行錯誤して， よーやくそれらしい estimator が得られた． わざわざここにその結果を書くほどでもないけど， あとから自分で参照するため， それから宮本さんなども不等分散な quasi-likelihood モデル使ってるので， ちょっとメモしておく．
• まず成長量 (でもなんでもいいんだけど) Gi の期待値を gi とする． この「ある個体は平均的にどれぐらい成長しそうか」 をあらわす gi は環境その他あれこれで決まるものとする． で， quasi-likelihood 推定のわくぐみでは， 成長量 Gi の分散 Var が gi の関数になってないといけない (たんなる定数でもよい)． ここでは gi の絶対値に分散が比例するとおく．

このように定めると個体ごとに定まる quasi-score，

を集団全体で最大化する 最大 quasi-likelihood 計算ができる (昨日しめしたように mlfitting や R の glm 使って)． この結果を使って σ2 を推定する． この推定計算においては最尤推定法を使えないので モーメント法で計算することになる．

…… でたぶんあってると思う． 「わりざんする」「平均値をとる」 というのに無意識の抵抗感があってなかなか思いつかなかった …… ということにしておくか． たぶん 1/N ではなく 1/(N-1) とすれば不偏だろう． こんなのを小標本系統でやるヒトがいるかどうかは疑問だが．
• 懸案の問題がかたづいたので 2020 研究室発． 北 12 生協で買いものなどして 2040 帰宅． 体重 73.6kg． 晩飯．
• 上の推定で間違ってないか気になったので， 数値例など作ってみてチェックを繰り返す． ま， 乱数標本で遊んでるだけですけど．
• [今日の素読]
  • あ， 今日は何もやってない
• [今日の運動]
  • こっちも
• [今日の食卓]
  • 朝 (0840): 米麦 0.5 合． ハクサイ・ニンジン・タマネギ・豆腐の味噌汁．
  • 昼 (1230): 低温研ちかくの北大北生協食堂で ブラック焼肉カレーなるあやしげな食物． 海草サラダつけて 514 円． 921 キロカロリー …… う， 若ものむけの食い物であったか． 蛇足ながらあまりおいしくなかった．
  • 晩 (2130): 米麦 0.8 合． ハクサイ・ニンジン・タマネギ・豆腐の味噌汁．

2002 年 12 月 15 日 (日)

• 0930 起床． 朝飯． コーヒー． でれでれ． 昼飯なども食ってから 1330 自宅発． 曇． 1340 研究室着．
• 研究室内のガスストーヴをつけると一酸化炭素で アタマがやられるのではと懸念して， つけずにいたところ …… やはり寒いとアタマが働かぬ， という原則はかなり普遍的かもしれんと 数時間後に確認できた．
• ということで記憶もとびがちなんだが …… 何やってたかな． Quasi-likelihood の文献を探したりしてみた． 北大・理・数学の図書室にある． 他に何をやったかというと …… それぐらいか? なんともはや．
• ということで何しに来たのかよくわからぬまま撤退． 2020 研究室発． 2030 帰宅． 体重 73.6kg．
• ひどく無為な休養日でした． 本屋とかにでも行けばよかった．
• [今日の素読]
  • うーむ ……
• [今日の運動]
  • ばて， てるわけでもないと思うが …… いやいや， ParseExcel & Quasi-likelihood ばて ということにしておくか
• [今日の食卓]
  • 朝 (1000): シリアル． コーヒー．
  • 昼 (1230): スパゲッティー． タマネギ・ピーマンの炒めもの．
  • 晩 (2130): 米麦 0.8 合． ハクサイ・タマネギ・ピーマン・油揚の炒めもの．

2002 年 12 月 16 日 (月)

• 0820 起床． 朝飯． 0850 自宅発． 曇． 0900 研究室着．
• なにかばたばたと計算してたんだけど思いだせん …… また quasi-likelihood 計算だったか?
• ひと区切りついたのでお茶部屋で休憩． M2 井田君と一般化線形モデルのギロン． なんか近ごろあちこちでこの計算方法 すすめてまわってるような気がする．
• かとー先生に MacOS X の挙動など教えてもらう．
• M1 大澤君と prunusdb のススめかたの相談． ともあれ現時点では Buzzaz 本にたたられていて， いそがしそうだ．
• D3 浦口さんのカエデ成長 quasi-likelihood モデルの ここまでの報告など． さーて， このハナシの統計モデリングは これで最後であってほしいところであるが ……
• 1250 北大構内走発． 積もった雪が踏み固められてしかもとけかけている． 気をつけないとすべりそうだ． 1325 研究室帰着． 昼飯． 昼飯後に大澤君のアメフト概論を拝聴する．
• 土曜日のぎょーむ日誌に TWINSPAN がどーのこーのと書いたら， もっとしっかり勉強せよとの御鞭撻と 関連情報を教えていただく …… ふーむ． この名義変数としか思えぬ優占度を 連続変数であるかのよーに扱うお作法， というのは ……
• 2000 年 3 月にメキシコのバハカリフォルニアで遭難した 中野繁さんの論文集「川と森の生態学」 (6000 円，ISBN4-8329-8021-1) がとどく． 中野論文を苫小牧の村上さんたち共著者たちが和訳してまとめたもの． 私もこれのお金だしてたんだなあ (すでに忘却)． 巻末用語集に関してコメントしたのは覚えてる． 出版社は 北海道大学図書刊行会 なんだけど， 同サイトにはまだ情報があがっていない．
• 1530 研究室出て理学部・数学科 (理 3 号館) 5F の図書室へ． かなり充実している． Quasi-likelihood 本 あった …… 書いた本人としてはすごく具体的に書いてるつもりなんだろうけど， 私からするとかなり数学よりな内容． ともあれ借りだしていく．
• ついでに農学部の図書室にも立ち寄る． ここの中央図書室はなかなか立派なんだが， 探してる Forest Science 古雑誌は無い． なんでも森林科学の図書室に収蔵とのことだが …… えーと北館 2F か． うすぐらく迷宮のごとき農学部内をさまよう． 森林科学図書室． はりがみ． 利用可能なのは火曜日・金曜日の 1200-1300 …… 週に 2 時間しか開いてないのか!! さすがは北大・農．
• 帰りがけに図書館本館にもたちよりごそごそと．
• もどってくると 1730 すぎ． 浦口さんから推定宿題メイル (本人は所用ですでに下校)． イタヤカエデはうまくいったけどオオモミジはうまくいかない? ふーむ． 一見するとオオモミジのほうがラクそうだけどな ……
• 推定計算やってみるとたしかにうまく収束しない． 謎． いろいろ試行錯誤してみる．
• その結果わかったことは， パラメーター推定の反復試行の途中において， 平均値がゼロ近くでぢたばたするとまずいことが生じうる， とゆー可能性だ． Quasi-score のマイナス側の極値が消滅して …… Newton-Raphson 法が破綻してるんだろうな．
• これを回避するには …… 数値計算上の都合として quasi-score 関数を某所でへし折ってしまうか． 一見すると乱暴なようだけど 今までやってることとは無矛盾であり， 疑似尤度を尤度そのものとみなして これを最大化するという方針とも合致している．
• さーて， 結果は …… この簡潔な処置で どうやらうまく推定できるようになったみたいだ． データのばらつきがでかい場合においても， いいかげんな初期値からスタートして 不安感なく収束するようになった． 一件落着．
• カエデの推定とは無関係なんだけど， かとー先生から R における ロバスト推定あれこれを教えていただく．
  • MASS ライブラリの の rlm 関数: これは M 統計量を使う方法のよーで
  • lqs ライブラリの の lqs 関数: これはどうも平均値のまわりの， 何と言うか「重みづけカーネル」 みたいなものを使って計算するみたいで …… しかし例をみるとうまく計算しているように見える
• 私はこの種の問題はぶーつすとらっぷ的に解決しよう， ともくろんでいたんだけど …… どちらが統計学的にはマシなんだろうか． 2110 研究室発． 雪どかどか． 2130 帰宅． 体重 73.6kg． 晩飯．
• [今日の運動]
  • 北大構内走 1250-1225． ストレッチング．
• [今日の食卓]
  • 朝 (0830): 米麦 0.6 合． ハクサイ・タマネギ・ピーマン・油揚の炒めもの．
  • 昼 (1335): 弁当． 研究室お茶部屋． 米麦 0.8 合． ハクサイ・タマネギ・ピーマン・油揚の炒めもの．
  • 晩 (2200): 米麦 0.8 合． ネギ・ニンジン・豆腐の味噌汁．

2002 年 12 月 17 日 (火)

• 0810 起床． 朝飯． コーヒー． 0900 自宅発． 雪． 除雪の結果として， 車道と歩道のあいだに「土手」が形成されつつある． 0910 研究室着．
• 先日おしえていただいた Ordination Methods for Ecologists なるペイジながめつつ計算のやりかたを検討してみる． なんとなく． やはりどういう統計モデルなのかよくわからん． まぁ， べつにすぐに問題点を整理しないといけないわけでもないし． ということで suspend．
• で prunusdb.pm のほうを resume して， ぢりぢりと取り組んでみる．
• カエデ quasi-likelihood 推定， じつはイタヤカエデのほうもうまくいってなかったようで． いやはや variance がでかいと収束計算しんどいのか? ともあれ， いまのところは昨日の修正 (?) quasi-score 法で うまく計算できてるようで．
• あ， 今日も弁当忘れてきた． 走るついでに食べに帰るか． 1230 北大構内走発． 1310 帰宅． シャワー． 着替え． 昼飯． 1350 研究室帰着．
• 研究室内の共用計算機のファイルを整理しろといった 指令がきてるな …… たしかになかなかひどい状況なのでなんとかするか．
• バックアップ方法など検討してみる． やはりネットワーク経由かな．
• とか調べていたはずなのに， いつのまにか甲山さんの G4 PowerBook の Eudora → AppleMail の変換作業をやっていたり. Eudora は腐っている．
• 1630 より Trendy セミナー ． 今日は低温研の加藤さんで八ヶ岳の Abies 稚樹のハナシ． 甲山さんの調査地のようにばたばた 林冠木が死ぬような場所ではないんだけど， 成長速度などは似ている． Abies PipeTree のモデリングでいままでよくわからなかった カサ型樹形個体の挙動とかわかった． これは PipeTree とは関係ないけど， 広葉樹の影響をどう評価するか， というあたりにまだ解析の余地ありだな．
• セミナー後に晩飯に出る． 今日から農学部で特別講義をやってる酒井聡樹さんも こちらに見えたのでいっしょに． たくさん食べてしまった．
• 2110 帰宅．
• [今日の素読]
  • なぜかススまぬ ……
• [今日の運動]
  • 雪の北大構内走 1230-1305． ストレッチング．
• [今日の食卓]
  • 朝 (0830): 米麦 0.8 合． タマネギ・モロヘイヤ・マイタケの炒めもの． ネギ・ニンジン・豆腐の味噌汁．
  • 昼 (1315): 弁当． 自宅． 米麦 0.8 合． タマネギ・モロヘイヤ・マイタケの炒めもの．
  • 晩 (1830): 北 15 の 中華美食楼 のコース (!) 料理． そう， いつもは定食しか食べないこの店で …… 棒々鶏． 双味蝦仁 (エビ二種類ソース)． 鮮蝦青瓜絡 (エビとキュウリせんぎり揚げ物 -- フライドキュウリうまい)． 紅焼鮑魚 (アワビとチンゲンサイ炒めもの)． 油林鶏 (からあげ)． 紅焼蹄筋 (スジ肉の炒めもの -- ねっとりした食感でおいしい)． 小籠包． 水餃子． 抜絲地瓜 (大学イモ)． この 8-9 人用コース 15000 円． みなさん青島ビールなど飲む．

2002 年 12 月 18 日 (水)

• 0810 起床． 朝飯． コーヒー． 0855 自宅発． 晴． 0905 研究室着．
• いろいろメイル書きなど． う． 苫小牧からの第二次攻撃予告．
• 今日は忘れものしてないように， という心づもりであったんだが Bazzaz 本を自宅においてきた． やれやれ． 北大構内走のついでに回収してくるか．
• いつものことながら， データ読みこみまわりを作るのは面倒だ． 下は個体の一部だけなんだが， これだけの項目があったりする．

	if ($head eq 'シュート数') {
		……(処理)……
	}
	elsif ($head eq '被摂食数') {
		……(処理)……
	}
	elsif ($head eq '葉巻き数') {
		……(処理)……
	}
	elsif ($head eq '葉枯れ数') {
		……(処理)……
	}
	elsif ($head eq 'マイナー') {
		……(処理)……
	}
	elsif ($head eq 'アリ来訪') {
		……(処理)……
	}
	elsif ($head eq '葉巻き') {
		……(処理)……
	}
	elsif ($head eq '葉枯れ') {
		……(処理)……
	}

まぁ， 大澤君のシウリザクラ観測データはやはりどちらかといえば， ふつーよりは複雑な構造をしてるということなんだろうなぁ．
• 1155 北大構内走発． 今日は早めに切り上げて 1225 帰宅． 昼飯など． 1245 自宅発． 1255 研究室もどる．
• 1300 より Bazzaz 本にも見所はあるはず再評価委員会 ． 今日の再評価委員は M1 大澤君． 他の人々が苦闘する Bazzaz ろぢっくを要領よくまとめている． さて，第 9 章 ``Ecological and genetic variation in early successional plants'' なんだが， 遺伝子解析はひとつしか登場しないのに， 遷移初期の草本 (結局この本はこればかり) 育成実験 (これは遺伝子なんぞ調べない) だけで環境不均質性と遺伝的多様性の関連について， まるまる一章ぶん論じきってしまうところは著者の面目躍如と言えよう．
• 1600 終了． そのあとお茶部屋で 「XP と Vine Linux 共存させてみてくださいよ」 という大澤挑発になぜか乗ってしまって Sharp Mebius なる ノート PC に 2 時間ほど費したのちに敗退してしまう． ここまで完全にヤラれてしまうのはひさしぶりだ． あとから調べてみたら (って先に調べろよ)， XP に不法占拠されている計算機上に Linux 人道支援を送り込むには 多少の手間と工夫 が必要とされるようだ．
• 撤退処理をしていると， 浦口さんが 「quasi-likelihood って AIC とかで モデル選択できないんじゃないですか」 という恐るべき追求を …… 言われてみたらそうかも． ああ， なぜ今まで気づかなかったのだろう．
• これまた今さらながらに調べてみると， quasi AIC (QAIC) なるものはある． しかしこれは下のような特殊な状況 (というか quasi likelihood 推定がもっとも使われる状況) にだけ対応するもののようだ．
  1. 一般化線形モデルで Poisson 分布 (離散・非負) の推定を行う
  2. 推定結果が overdispersion (平均 < 分散) であれば quasi likelihood 法で 推定しなおす
  3. 推定を 2. のように行った場合には (overdispersion を勘案した) QAIC でモデルを評価する (quasi likelihood 推定でないときは普通の AIC)
• しかしいまやってるのは 「指数関数族のひとつの分布であって， 平均と分散の関数関係だけがわかっている」 という状況のもとでの推定なんだよね． うーむ． 擬似尤度そのものはパラメーター推定には使えても， しょせんはニセ尤度なんで確率論的モデル選択できん， ということのようで …… なるほど (Poisson 分布問題以外では) 擬似尤度が普及しないはずだ．
• もうちょっと調べてみると …… さすがわ統計学者ども． なにやら quasi-likelihood で モデル選択ごとき操作をやろうとぢたばたしてるようだ． しかしこのあたりは理論的にも応用的にも まだ確立してない分野のような気がする．
• Quasi 文献はまた図書館にいってみることにして， 別の策を考えてみる．
• …… うーむ， ひょっとすると， 以前にほうりだしていた構想を使えるのではなかろーか． それはこういうややこしい確率論的モデルで ……

  (観測された値) = (ある確率分布の標本) + (測定誤差)

  右辺の第一項は自然がつくり出した部分 (これを推定したい) で 第二項が人間によってつくり出された部分．
• で， これを満たすような尤度方程式を書いてみる． やはり面倒というか時間のかかりそうな積分が入るか． しかし他に手がなさそう．
• これを mlfitting に計算させる plugin_function.cc 書く． しかし試験運転などやったりすると眠れなくなりそうなんで， コンパイルが通ったところで撤退． 2420 研究室発． 2430 帰宅．
• ここ二日間ほどばたばたしていて， あちこちにメイル返事が書けづ． とどこおっておりますです．
• [今日の運動]
  • 北大構内走 1125-1225． ストレッチング．
• [今日の食卓]
  • 朝 (0840): 米麦 0.5 合． タマネギ・モロヘイヤ・マイタケの炒めもの． ネギ・ニンジン・豆腐の味噌汁．
  • 昼 (1230): 米麦 0.5 合． タマネギ・モロヘイヤ・マイタケの炒めもの． ネギ・ニンジン・豆腐の味噌汁．
  • 晩 (2120): 弁当． 研究室お茶部屋． 米麦 0.8 合． タマネギ・モロヘイヤ・マイタケの炒めもの． ネギ・ニンジン・豆腐の味噌汁．

2002 年 12 月 19 日 (木)

• 0820 起床． ねむい． 朝飯． コーヒー． 0855 自宅発． 晴． 0905 研究室着．
• いろいろ進行中の計算問題のうち， 何をやったらよいのかわからなくなったので， 「ぎょーむ日誌」を書いてみる． ようやく昨日の作業進捗を思い出してきた． しかし日誌書きは時間かかる．
• とりあえず北大図書室めぐり． まず経済学部に行って generalized linear model の教科書を二冊借用する． うち一冊は標準的な教科書である McCullagh & Nelder, 1989． 次に数学科に行って先日の quasi likelihood 本を返却． なぜか文学部に Biometrika の少し前の巻がそろっているので， そちらに潜入して quasi deviance criterion 論文の複写．
• 複合確率分布の最尤推定問題のつづきにとりくむ． まず浦口さんに教えていただいた Math::Random (これは CPAN に含まれていて，かなりの種類の乱数を生成できる) を使って， Gamma 分布 + Gauss 分布な標本集団を生成してみる．
• で， これを最尤推定プログラム mlfitting と複合分布用 plugin_function.so に与えてみると …… 試験運転してみると予想どーりうまくいかない． さて， ということで問題が発生してる箇所を特定していく．
• なんだ数値積分の最初の最初でコケてるのか， とわかったところで時間切れ． 昼飯．
• 1300 より研究室の歳末掃除． 私はお茶部屋掃除に配属． 床をみがきゴミをどんどんだしていく．
• 1600 すぎに終了． またなぜか甲山さんの Eudora → AppleMail 移行問題に取り組む． MacOS X の「アドレスブック」とやらが データ入出力にもちいる vCard なる書式はわかったので， Perl で変換スクリプトを組む． おお， MacOS X で vi とか当り前のように使えて まるっきり Unix 系 OS だ …… しかし大半の MacOS X ユーザーはかかる便利ツール群の存在 にすら気づかずに生涯をまっとうするのだろう．
• で， 書式変換はうまくいった． しかし日本語文字がばける． このアドレスブックなるソフトウェアは いったいどういう文字コードなら正しく読解してくれるのか? ShiftJIS だめ JIS だめ EUC だめ Unicode だめ …… なぞだ．
• やはりこちらも時間切れ． 1730 講座のみなさんうちそろって忘年会に出発． 1745 JR 札幌駅東側の うおや一丁 に到着． 忘年会．

[忘年会幹事長] 北日本における忘年会の幹事正装．

• 2000 にここは終わって， 二次会はなんとなくとなりの居酒屋たじま屋に． 損害報告概要 …… ひどく酩酊してもの投げたり隣人の首をしめてまわったりした D1 一名． 行方不明になったんで捜索してみると トイレで沈没してる状況で発見された M2 一名． 「あと 5 年間結婚しません」 と力強く宣言してから 前後不覚に泥酔して昏睡してしまった女性大学院生一名．
• 2340 同発． 大学にもどる人数が多そうなので， 暴れまわる酩酊院生の輸送はそちらにおまかせしてそのまま帰る． 2400 帰宅．
• [今日の食卓]
  • 朝 (0840): 米麦 0.5 合． ネギ・ニンジン・豆腐の味噌汁．
  • 昼 (1230): 研究室お茶部屋． パン． コーヒー．
  • 晩 (1800): 忘年会であれこれと．

2002 年 12 月 20 日 (金)

• 0820 起床． コーヒー． 0855 自宅発． 曇． 0905 研究室着． 嵐のあとの静けさ．
• 積分やりながら多次元微分方程式を解く 複合確率分布の最尤推定問題の続き． これで午前中おわった． ぢりぢりとしかススまず．
• 昼飯． 今日はなんとなく運動休養日になりそう．
• なぜかまた甲山さん G4 PowerBook MacOS X 上の Eudora → AppleMail 移行問題． いつものごとくかとー先生がナゾを解明した． 呪われ vCard 書式の 1 アイテムは例えばこう書けばよい．

	BEGIN:VCARD
	FN;CHARSET=UTF-8:久保拓弥
	EMAIL;WORK:kubo@ees.hokudai.ac.jp
	END:VCARD

CHARSET で明示的に指定すればどんな文字コードでも OK なのかもしれない (しかしそういう実験はやってない)． で，変換の手順はこうなる．
1. 腐れ Eudora から何としても nickname というファイルをテキスト化する (これが Eudora ヴァージョン依存なので 「こうせよ」と断定できない)
2. MacOS X の「Unix 側」はまったくもって Unix なので 上で変換したテキストファイルの 行末文字を一括置換しないといけない …… tr '\r' '\n' (しかしこれぐらいのことは次の Perl スクリプト内にて処理すべきだったのではなかろうか)
3. 全データを Perl スクリプトで上のような BEGIN:VCARD …… END:VCARD 書式に一括変換する．
4. はなはだ腹だたしくも， ここで MacOS X のテキストエディター による手動作業 …… 文字コードを UTF-8 に変換する (あとでかとー先生から Unix の iconv 使えばいいのにとの御指導をうける)
5. 呪われ vCard の行末文字はなぜか M$DO$ 形式でないと いけないんでまた一括置換 …… sed 's/\n/\r\n/g' (これも Perl 内でやるべきだよなぁ ……ひとつ前でよけいな操作してるんで， そうはできなかった)
やれやれ， というほかない作業に決着．
• M1 大澤君の Sharp Mebius における XP と Vine Linux 共存問題． 私が一昨日に boot 不可能にしてやった XP 領域を すばやく復旧してきたんで …… USB FDD を準備して Vine Linux インストールに再びとりくむ． 要点はただひとつ， booter (ここでは LILO) を MBR に 書き込んではいけない ということに尽きる．
• で， 先日みつけた 共存ペイジ を参考にしながら面倒な作業をやろうとしたんだけど …… この Sharp Mebius に関しては， 結果としてここに書かれてるような面倒なコトはやらなくてもすんだ．
• HDD は一個だけで，先頭からこういうふうに切られているわけだが ……

	MBR        WindowsXP に支配されてる
	/dev/hda1  WindowsXP       (16 GB)
	/dev/hda2  Linux swap 領域 (0.2 GB)
	/dev/hda3  Linux           (3.8 GB)

Vine Linux インストール時に LILO は LBA32 強制指定して かつ /dev/hda3 の先頭に書き込ませた． こんなところに booter を書き込んでも機械のほうは自力では 見つけてくれないはずなんだけど (ということで，boot FD あるいは XP booter に起動させる 予定だった) …… なぜか Mebius の BIOS が /dev/hda3 アタマの LILO を勝手に見つけてくれるようになったんで (LBA32 強制指定と関係アリ?)， これでつつがなく Linux 人道支援を達成できた． めでたしめでたし．
• というようなコトは片づいたわけだが， 数値的最尤推定法問題は進展せず． 撤退． 1840 研究室発． 1900 帰宅． 体重 73.4kg． 忘年会太りをまぬがれているんだろうか．
• 晩飯食ったあとに， 過去三日間ほど対応できなかったメイルに まとめてお返事書き．
• しかし計算仕事はススまず．
• [今日の素読]
  • こっちもススまず．
• [今日の運動]
  • まぬけな運動休養日．
• [今日の食卓]
  • 朝 (0840): 昨晩の忘年会で食べすぎて胃のあたりがキモチ悪いんで コーヒーだけ．
  • 昼 (1250): 研究室お茶部屋． パン． コーヒー．
  • 晩 (2010): 米麦 0.8 合． ハクサイ・ネギ・シイタケ・卵の炒めもの．