ぎょーむ日誌 2002-12-13

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 合． ハクサイ・ニンジン・タマネギ・豆腐の味噌汁．