ぎょーむ日誌 2005-02-(01-10)
2005 年 02 月 01 日 (火)
-
0710 起床．
朝飯．
コーヒー．
0850 自宅発．
晴．
0900 研究室着．
-
D 論発表会二日目．
今日は地域計画のほうのバハさんから．
家の値段の spatial autoregressive model．
これは先日
しらべた Ord (1975) の方法で最尤推定できる問題になってた
(後記: じつはこれが解けるかどうかは
計算してみないことにはわからない，
と翌日になって理解できた)
……
まあ，
あちらの分野の定番ワザなのだろう．
しかし標本数 143 なのにパラメーター数あまりにも多すぎるよーな
気がしたのでそのあたり質問．
1010 終了．
生態のほうは午後からである．
-
昼飯くってから午後の発表に出席
……
浦口さん・平尾君が無事に発表おえて 1410 おひらき．
えーい，
しかし本日はいつもにもましてまぬけ質問などやってしまった．
反省．
とくに平尾君の集団遺伝学ハナシ聞いてると，
アタマの中で M2 立花さんのミミコウモリ集団遺伝学の難問どもが
励起されてしまって，
ですね．
-
で，
D 論発表聞きにこられていた森林総研の永光さんと
お茶部屋でそのあたりの相談を．
すると，
当事者たる立花さんから
「ゑくせる上の手動計算まちがい」
の報告が
……
またしてもゑくせる禍ですよ．
-
で，
計算結果のかわったしょっくも乗りこえて，
いろいろ 3 人で議論してみたんだが
……
うーむ，
というかんぢ．
-
D 論発表会聞きにこられていた鍋嶋さんと苫小牧直径成長モデル相談少々．
帰無仮説的に年変動をランダム効果にすることはできる
……
計算かえって面倒になる，
というのはともかく，
この二重に混合なモデルが現象をよく説明できてしまいそうな気がして，
なかなかオソろしい．
-
森林総研北海道支所に短期滞在研究に来ている
下野さんに筑波大学の様子などうかがう．
-
北大ネット管理者から全ユーザーにばらまかれたまぬけメイル．
かれらはハードディスク増設するだけのために，
全ユーザーのパスワードを「初期ぱすわーど」
に戻さねばならんらしい．
どこのど素人に管理者やらせてるんだか．
いやはや．
``in Japanes'' も原文ママ．
------------------------ in Japanes ----------------------------
HINESメール利用者 各位
メールサーバ移行に伴う初期パスワードについて
日頃は、HINESの運用に対しご協力いただき感謝申し上げます。
さて、現在使用しているメールサーバを新しいメールサーバへ
移行（16年度中、個人容量の増量等予定）するべく準備をしています。
それに伴い、新メールサーバでは、PID・転送設定等に変更はありま
せんが、システムの都合上パスワードを初期パスワードにて設定する
ことになりました。
そこで、初期パスワードが分からない方は、このメールにてPID(*1)
・所属(*2)・内線等電話を必ず明記の上、2月13日(日)（以降も受け付け
ますが）までに返信してください。PID・初期パスワードの記載された
用紙を学内便等にて本人宛へ送付いたします。
なお、移行の詳細につきましては、決まり次第メール・HINESのホーム
ページ等にてお知らせいたしますので、もうしばらくお待ちください。
よろしくお願いいたします。
注意：*1・PIDが分からない場合は、このメールのヘッダー部のToに
設定(To: xx123@pop.・・・ のxx123がPID)されていますので、
ご確認ください。
*2・名誉教授の方について所属は、北大時の部局名と郵送先が
部局以外の場合は、郵送先住所をご記入ください。
・返信は、初期パスワードが分からない方だけお願いします。
追伸：メールサーバ移行に伴い、来年度への継続作業ができませんので、
今年度は、継続の手続きを行わない予定です。
(英文らしきもの省略)
2005.2.1
<<<<<<<<<<<<<<<<<<<<<<<<<<<<
北海道大学情報基盤センター
ネットワーク係
011-706-3456/2944
hines-info@iic.hokudai.ac.jp
>>>>>>>>>>>>>>>>>>>>>>>>>>>>
で，
さっそく「自分の初期ぱすわーど教えてください」
と北大全体にメイルを送信してしまうやつがあらわれるし
……
1910 研究室発．
1930 帰宅．
体重 74.4kg．
運動．
晩飯．
う．
混合モデルに関する難しそうな質問メイルが
……
[今日の運動]
[今日の食卓]
- 朝 (0730):
米麦 0.7 合．
ハクサイ・ネギ・ナメコ・ワカメ・コンニャク・豆腐・煮干の味噌汁．
- 昼 (1150):
研究室お茶部屋．
米麦 0.7 合．
ハクサイ・ネギ・ナメコ・ワカメ・コンニャク・豆腐・煮干の味噌汁．
- 晩 (2230):
米麦 1.0 合．
ハクサイ・ネギ・ナメコ・ワカメ・コンニャク・豆腐・煮干の味噌汁．
エノキダケ炒めもの．
2005 年 02 月 02 日 (水)
-
0800 起床．
朝飯．
コーヒー．
0910 自宅発．
晴．
0920 研究室着．
-
いろいろなヒトたちのスケジュールのかねあいから，
ずいぶんと早い時期に開催される「追い出しコンパ」情報めも．
ここに提示する意味はないんだけど，
もしかしたらこういうところに置いといたほうが便利かも，
という実験
……
というほどでもないか．
[日にち] 2月16日（水）
[会場と時間] 一次会 「ひだまり庭」 19:00-21:00
住所：北11西3
tel：011-709-5850
二次会 「がらり」
住所：北18西4
tel：011-707-7155
※地図をお茶部屋に掲示します。
[会費] 4000円（一次会）
2300円（二次会）
ここしばらくの修論・D 論発表など見ていると，
R
の使用とともに，
いよいよ一般化線形モデル (glm) も普及してきた．
その「使いたおし」ぐあいはヒトそれぞれで，
単に「なんでも『分散分析』」が「なんでも glm()」
になっただけのもいるし，
その一方で
「統計モデルは現象を表現する手段である」
という理解にもとづいて
適切なモデルを構築できる大学院生たちもいる．
そのなかで特に賢明なる大学院生たちは，
単純な一般化線形モデル (GLM)では
十分に説明しきれない現象があることに気づき
(個体差問題など)，
そのあたりを一般化線形混合モデル
(GLMM) で何とかできるかも，
と模索している
(ただし混合モデル化が常にこのあたりの解決策というわけではない)．
一年ほど前に R で一般化線形混合モデル
のパラメーターを推定計算する関数として
glmmML() を紹介
してみたんだけど，
このあたりの現状について簡単にまとめてみよう．
まずは最尤法にもとづく推定用の関数:
-
glmmML()
(in glmmML package):
Göran Broström 作．
CRAN に入ってる．
推定結果に
stepAIC()
を適用できない
(で，
私は自作のモデル選択関数でそれに対応している)．
-
glmm()
(in repeated package):
Jim Lindsey 作．
CRAN に入ってない．
推定結果に
stepAIC()
を適用できる．
AIC なども自動的に計算してくれる．
後者は苫小牧の平尾君に教えてもらった．
最尤法を使う
利点は本モノの尤度を使っているのでモデル選択などができる，
ということ．
欠点は，
現時点ではランダム効果をひとつしかいれることができない
(しかも定数のみ)，
推定計算の収束に失敗する可能性が PQL 法より高いかも，
というあたりか．
次に罰則つき擬似尤度 (penalized quasi-likelihood; PQL)
法にもとづく推定関数 (どちらも CRAN に入っている):
-
GLMM()
(in lme4 pacakge):
Douglas Bates
and Deepayan Sarkar 作．
-
glmmPQL()
(in MASS pacakge, nlme package も利用)
これらの相違はふたつある: 一つはオブジェクトの属するクラスで
S4 class (GLMM)
か S3 class (glmmPQL) か
……
まあ，
これは一般ユーザーにはあまり関係ない;
第二点は計算結果の一部が両者で異なる
(計算方法が異なるため)．
PQL 法を使う利点は簡単に multi-level モデルの推定計算ができる，
というもの．
欠点は PQL は尤度ではないので，
尤度にもとづくいろいろな方法 (モデル選択など)
が使えない，
というところである．
GLMM()
や
glmmPQL()
を使った計算結果オブジェクトには最大化対数尤度や AIC
(後者のみ) も含まれているけれど，
これらは最尤法で得られる尤度とは別モノである．
これら PQL な GLMM 計算関数は，
「一般化」ではない線形混合モデル (LMM?)
を解く関数を利用して擬似尤度を最大化している．
線形混合モデルというのは，
例によって何もかもが等分散正規分布の世界である．
しかしこの仮定をうまく利用することで
多重階層なモデルの計算が可能になっている．
最尤法を使って multi-level なモデルの推定やりたければ
……
自分で関数をかくしかないと思う．
この場合，
重積分はしんどいので，
MCMC 法などモンテカルロ法を援用したほうが良いだろう．
あるいは，
どうせ MCMC 使うならもはやベイズ推定で，
というふうになってくるのかもしれない．
1000 フィールド科学センターの
前川さん
の院生の渡辺君の統計学こんさる．
私が
先日
の修論発表会であれこれと統計学的難クセをつけたため．
データの基本構造について教えてもらい，
魚セット (メス 1, 回遊型オス 1, 残留型オス 1)
の「個性」を考慮した一般化線形混合モデル (GLMM)
で何とかなりそう，
という方針を検討する．
今日は生データなく
(生態学の悪しき伝統のひとつというべき)
「わりざん値」
な情報欠落ずみデータだけだったんで，
また後日に解析をススめることに．
北大生協から
大阪大会
の航空券と内地入国の査証が届いたので，
カードで支払う．
領収書は三浦さんにお渡しして処理をお願いする．
昼飯．
なんとなく統計学勉強モードに．
GLMM など混合モデルでのモデル選択規準は，
単純にランダム効果のパラメーター数くわえて良さそう，
というあたりの再確認
……
まあ，
Akaike's Bayesian Information Criterion
(ABIC) なんかがそういうカタチになっているから，
というのだけでは正当化もへちまもないわけだが．
さらに統計学勉強が続いてしまう．
昨日書いたことにまたたまた一部まちがいあったわけで
……
バハさんの simultaneous autoregressive (SAR) モデルの推定計算，
あれで最尤推定になってると思ってたんだけど，
そうではなく pseudo likelihood estimator にすぎなかった．
難しい．
ある限定された SAR で最尤推定を行なう
Ord (1975) の方法というのを原論文
Ord, K. 1975.
Estimation methods for models of spatial interaction.
Journal of the American Statistical Association 70: 120-126.
を読みなおして
(以前はいいかげんに眺めていただけだった)
勉強してみる．
空間相関のパラメーターを推定するのは，
それほど簡単ではない．
ひとことで言えば，
距離荷重行列の固有値をすべて計算して
……
となるわけで，
複素固有値がまじってきたらどうすんだと思うわけだが，
Ord 先生いわく絶対値が大きいところだけ取ればよい，
とのこと．
まあ，
絶対値の大きい固有値が複素数になったらこの方法はそれでもうダメ
になるわけだが．
小林さんがアカマツデータ解析をしつつ，
論文の準備にとりくんでおられるようなので，
とりあえず依頼された LaTeX ファイル変換して先方に送ってみる．
簡単な Perl スクリプトでフィルタリングしただけのものなんだけど．
で，
けっきょく Ord のアルゴリズムとは
-
SAR モデルを空間依存項
(推定すべき空間相関なパラメーター
ρ
を含む)
とそうでないところに分割する
-
空間に依存してないところ
Y = X β + u
について
ordinary least square (OLS) を適用して
「残差」
u
を求める
-
その
u
を使って
u = ρ W u + ε
から
ρ
の推定計算するわけだが，
このときに固有値計算が必要になる
-
3. の結果にもとづいて
2. に必要な新しい
Y
を計算する
-
値が変化しなくなるまで 2. -- 4. をくりかえす
というもののよーだ．
まあ，
ふつーのユーザーはこういう問題に出あったら
(以前に苫小牧の平尾君にご教示していただいたとーり)
R
の
spdep
package の
lagsarlm()
関数でもって計算してしまえばよいわけだが．
これはまさに上述のアルゴリズムを実装した関数である．
1920 研究室発．
1930 帰宅．
晩飯．
明日の輪読会の予習．
いやはや，
なかなか手ごわい内容なんで．
上記「追いコン」幹事長として東奔西走する尊敬すべき
M1 女性大学院生殿から教官どもあての
「嘆願書」と血書したメイルいただき，
「慢性的財政難の学生だけでは、賄い切れません」
とのこととて
「二人分の料金を頂きませんことを、
此処に、切にお願い申し上げる次第でございます」
なる清廉潔白というほかない直訴である．
空腹のあまり勉学すらままならぬ
大学院生たちのまさにまさに窮乏のきわみというべき現状に袖を濡らしつつ，
御心配には及びますまいとの返信したためるうちにも，
以下のような見事なるしめくくりの
感銘のあまりの引用についてのお許しいただかずにはいられないのであった．
なお、当日は出発前に、集金に上がらせて頂く所存でございます。
快く承諾して下さる皆様の懐の深さを信じて疑いません。
[今日の運動]
-
腹筋運動 30 ×
3 回．
腕立ふせ 5 ×
3 回．
[今日の食卓]
- 朝 (0830):
米麦 0.5 合．
ピーマン卵炒飯．
- 昼 (1320):
研究室お茶部屋．
米麦 0.7 合．
レトルトパウチドカレー．
- 晩 (2130):
米麦 1.0 合．
タマネギ・マイタケ・豆腐・カツオブシの味噌汁．
モヤシ・ニラの炒めもの．
2005 年 02 月 03 日 (木)
-
0730 起床．
朝飯．
コーヒー．
0820 自宅発．
晴．
0830 研究室着．
輪読会の予習．
-
輪読会中止．
えー，
「またしても」というか．
今回は風邪だそうで．
Life History Evolution in Plants
，
「さーて，今日は第 8 章か」
なる状況はすでに三度目なんだけど，
そのたびに交通事故とか親族に不幸あったりとか風邪をひいたりとか
……
なンか呪われてますよねえ．
-
時間あいたので
ぱいぷ樹木原稿みなおし．
窓の外は雪．
-
昨日の「嘆願書」について議論が生じている．
文言のなかの
「二人分の料金を頂きませんことを、
此処に、切にお願い申し上げる次第でございます」，
この「頂きません」とは何なのだろうか?
という解釈問題で皆さん苦しんでおられる．
「嘆願書」全体および前後の文脈から，
このあたりわかりやすく意訳すれば
「つべこべいわず二人分しはらいやがれ」，
つまり幹事長殿からの
厳粛なる指示であることは疑いの余地ない．
しかしながら，
なぜそれが「頂きません」
と表現されるのか，
その文法的な側面についての検討も必要だろう．
-
これって北海道弁では，
という説もあったんだけど，
それだけではないのかもしれない．
近い表現を考えてみると，
これは「頂戴せん」だと考えればつじつまがあうような気がする．
「ん」は主体的意志をあらわす「む」の音便で，
意味としては
「頂戴してやろう」
となる．
おそらく，
北海道の一部では
「『頂戴』す」を「『いただきま』す」と読むのだろう．
「(私が) 二人分の料金を頂戴せんことを」
とすれば (これだけ見れば) 意味がつうじそうだ．
しかし，
そうだとすると
「お願い申し上げる」に続けるのはヘンだ．
あるいは，
このあたりに研究室の実質的な権力者である女性大学院生たちの
政治的巧妙さがあるのかもしれない．
-
東アジアにおける権力のありようとは，
被支配者たちが自発的に王に従うことをもって最上とするものだろう．
昨日の引用文にはっきりと示されているように，
支払う側であるわれわれには選択の自由など
そもそも全くないのだけれど，
そうであっても，
いやそういう場合にこそ「自発的に支払う」
というカタチを取らねばなるまい．
「お願い申し上げる」うんぬんはまさにそのための修辞だ．
では，
「二人分の料金を支払っていただけますよう
……
お願い申し上げる
……」
というふうに書かないのはなぜか?
-
これはおそらく「特殊な読み手たちへの配慮」だろう．
上のように丁寧に書いてしまうと，
文言を「字句どおり」にしか読めない理系男どもは
「あ，自分には選択の余地あるんだ」
と錯覚してしまうかもしれないではないか．
それを
「二人分の料金を頂きませんことを」
と，
あえて，
わざと発話者の主体的意志を (北海道的に) あらわす
「頂きません」
を変則的にもちいることで支払い側の利己的な誤読をふせごうとした，
と思惟される
(おかげで私は文意を正しくとれた)．
企図は明確であったものの，
手法がかくも複雑であったために一部に混乱が生じてしまった．
あるいは民草たちが権力者の意志を推し量りかねて右往左往する
この状態こそ，
その真のねらいだったのかもしれない．
-
……
というようなことを検討議論しつつ昼飯．
書いた本人の談話
「古文調で書いてみました，えへへ」．
-
月末にまた試験監督．
二次試験．
-
甲山さんがまたぱいぷ樹木原稿のなおし提案もってきてくれたんで，
それを反映してみる．
もう一回よみなおしてみる
……
1700 までかかって作業終了．
その後の処理あれこれ．
-
1830 研究室発．
2000 帰宅．
運動．
晩飯．
-
[今日の運動]
-
エアロバイク 50 分間．
-
腹筋運動 30 ×
3 回．
腕立ふせ 5 ×
3 回．
-
[今日の食卓]
- 朝 (0750):
米麦 0.7 合．
タマネギ・マイタケ・豆腐・カツオブシの味噌汁．
コマツナ．
- 昼 (1320):
研究室お茶部屋．
米麦 0.7 合．
タマネギ・マイタケ・豆腐・カツオブシの味噌汁．
- 晩 (2220):
米麦 1.0 合．
モヤシ・ニラ卵炒飯．
コマツナ．
2005 年 02 月 04 日 (金)
-
0700 起床．
朝飯．
コーヒー．
0840 自宅発．
曇．
0850 研究室着．
-
ネット雑用．
進化植物研究会 (今年も夕張で開催)
のペイジ，
もうすぐ公開の予定．
-
う．
R
上で OpenGL な三次元作図のできる
rgl
で遊んでしまった．
こりゃ，
便利だわ
……
-
しばらく放置してたアカマツ論文原稿をみなおす．
現在，
小林さんのほうで光合成に関するほうの論文をまとめる準備をされているので，
当方はその残骸というべきアヤしげシミュレイションについて，
廃品再利用の可能性を検討してみる．
とりあえず，
-
光合成まわりは小林さん論文で使うので，
こちらからは外す
(じつはシミュレイションと光合成測定値のあいだに
何の関係もない，
というのがこれまでの弱点のひとつだった)
-
新しい投稿先にあわせて内容の改定
うーむ，
観測データは利用してるものの，
ひどくモデルな方向に偏った内容なんで
……
なんとか
Ecological Modelling
あたりに押しこみたいところではあるんだが．
とりあえず，
Elsevier のサイトなど調べる．
LaTeX file guidelines
あるのはよろしい．
-
昼飯．
昨晩，
大学院生たちがおでんパーティーやったんで，
その残りをいただく
……
というか，
これからしばらくおでんばかり食い続けることになるかも．
-
甲山さんとインドネシアデータ解析相談．
なンかよくわからんデータなんで，
てきとうにガンマ分布あてはめることに．
-
農学部の
造林
研究室から城田さんがとつぜんみえる．
しべりあなカラマツ林の開空度計算な相談．
まあ，
先日考えたよーな高速化ライブラリも必要ではあるが，
従来型の計算ライブラリも保守せんといけませんな．
-
お茶部屋で他のヒトたちもまじえて雑談．
-
露崎さんが
洞爺湖の中島
に生じたナゾの穴のなかでエゾシカ食害をまぬがれてる
シダのデータを調べ直しておられるよーで．
-
またアカマツ原稿のみなおし．
現時点の Introduction は二部構成とでもいうべき長大なモノで，
小林さんが書いた汚染物質曝露下の樹木に関する前半と，
私がでっちあげたそのあたりのシミュレイションに関する
先行研究の概要あれこれである．
これ，
シミュレイションだけの論文にするなら，
前半を書き改めないといかんな．
そもそも，
この前半部は小林さんが光合成論文のほうで使うだろうし．
-
1900 研究室発．
1910 帰宅．
運動．
晩飯．
-
[今日の運動]
-
[今日の食卓]
- 朝 (0730):
米麦 0.7 合．
トマト・豆腐のスープ．
- 昼 (1300):
研究室お茶部屋．
米麦 0.7 合．
おでん．
- 晩 (2110):
米麦 1.0 合．
おでん．
これは研究室から持ちかえったもの．
2005 年 02 月 05 日 (土)
-
0820 起床．
外はかなり雪がふってる．
朝飯．
コーヒー．
洗濯．
提出用ぱいぷ樹木原稿関連ファイルなあれこれ処理．
-
1040 自宅発．
雪．
1100 研究室着．
怠業．
昼飯．
北大生協書籍部まで雪中散歩．
-
統計学勉強．
多変量正規分布 (multivariate normal distribution)
の条件つき分布について．
他の変量が決まってるという条件のもとで，
ある変量の確率分布は
(共分散行列で決定される)
ある定数分散の正規分布となる．
うーむ
……
これは私のへっぽこというかアテにならぬ思い込み
(分散も増減すると思ってた
--- 平均は当然ながら与えられた条件によって増減するわけだが)
がまた間違ってたわけで．
いやはや．
-
蛇足ながら
R
で逆行列を計算する関数は
solve()
(これはちょっといやはやな名前ですな)．
あるいは MASS package には ginv()
なる
ムーア・ペンローズ型一般化逆行列を使った
計算関数もある．
-
ということで，
Markov Chain Monte Carlo (MCMC) 法において，
多変量正規乱数の Gibbs sampler はとても簡単で
計算負荷の少ないものとなる，
ということで．
-
そして，
空間相関ある GLMM みたいな状況で使う場合には，
線形予測子全体が多変量正規分布になる
……
と考えたほうが，
やはり私には理解しやすい．
理解しやすい，
というかそういうふうに計算しないと，
ダメだ．
Random effects の平均がゼロにならない，
というのはなんだかヘンなかんぢではあるけれど．
いや，
集団全体ではゼロになってんのかな．
-
ともあれ，
このような技法でもって苫小牧寄生蜂問題だの，
林冠動態ギャップ拡大問題
(こちらは私にとっては 10 年来の難問だった)
だのといった愉快な面倒は解明されていくのでありましょう
……
とゆーよーなハナシを「空き時間の埋め草」
として
進化植物研究会
(つい先ほど手直しして公開)
でやらんといかんかな，
と思っていたんだけど，
西村さんの尽力で大原研のかたがしゃべってくれそうな
(現時点ではまだ演題などいれてない)
……
ということで，
道内の生態学研究者たちが
MCMC 洗礼をあびるのは先のばしされた．
もとをただせば小菅君が小菅人脈
(これは旭川の陸自第二師団内部から
故郷たる横須賀の近所の総理実家その他にまでつながってるらしい)
を
「こんなことのために貴重な人脈を使うのは」
と出し惜しみしたのが，
今回の混乱のそもそもの発端といってよいだろう．
そうにちがいない．
-
1900 研究室発．
雪．
1910 帰宅．
晩飯．
早めに寝てしまったものの，
夜中に目覚めてしまった．
あいかわらずの降雪．
-
[今日の運動]
-
[今日の食卓]
- 朝 (0850):
米麦 0.5 合．
トマト・豆腐のスープ．
- 昼 (1230):
研究室お茶部屋．
米麦 0.8 合．
おでんの残り
……
には固形物が含まれてはいなかったのだが．
- 晩 (2000):
米麦 1.0 合．
トマト・豆腐のスープ．
2005 年 02 月 06 日 (日)
-
0820 起床．
朝飯．
コーヒー．
-
ぱいぷ樹木
のソースコードを少しいぢる．
変数名をパイプ樹木原稿のそれにあわせる作業．
Perl スクリプトを使って一括置換．
これは何の問題もなくできたんだが
……
動作が今までのと異なる．
謎．
あちこち調べてみる．
-
ちょっと解明できんので，
これはまたあとで登校後に．
北大構内走にでようとしたものの，
視界とざすような降雪．
しかたなく室内運動．
体重 73.8kg．
うーむ．
洗濯．
昼飯．
-
1600 自宅発．
雪．
1610 研究室着．
-
ここしばらく
ssh
で remote login するときにパスワード要求なんぞにナヤまされてたんだけど
($HOME/.ssh/authorized_keys その他はきちんと設定してる)
……
ふと思いついて server 側の
/var/log/secure
のぞいてみると
Authentication refused: bad ownership or modes for directory /home/kubo
なる sshd 殿からの怒りの message がずらずらずらと並んでる．
はて，
面妖なといぶかしみつつ，
chmod 711 /home/kubo
で解決．
chmod 700 /home/kubo
でもよい
……
が，
こちらでは不便な場合もある．
-
ぱいぷ樹木シミュレイター
PipeTree
の挙動不審問題，
表面的な原因だけは特定できた．
ぢつにおどろくべきものだった．
おそらく，
コンパイラー gcc の version 依存である．
新コードを gcc-c++-2.95.3-2vl22
(in Vine Linux 2.6r)
でコンパイルすれば旧来どーりの動作であり，
いっぽう，
旧コードを gcc-c++-3.3.2-0vl8
(in Vine Linux 3.1)
でコンパイルすると，
やはり午前中に観察されたのと同じおかしな動作になる
……
くそう，
呪われ言語め．
-
おかしな動作とは，
具体的には樹齢 3 ぐらいで梢端がスッとんで個体が死んでしまう，
という奇怪なしろものである
(梢端死が個体の死，
というのは私が設定したことなんだけど，
どうしていきなりこの leader shoot が突然死するのか
その理由がわからない)．
こんなことが，
コンパイラ (あるいは可能性ひくいものの動的ライブラリーの何か)
に依存して発生するとは．
現時点の憶測としては，
standard template library (STL) か名前空間まわりだろうなぁ，
という気がする．
-
コンパイルそのものは問題なくできる
……
なぜかというと，
面倒になりそうな部分については，
1.5 年ほど前
に解決ずみだったんで．
よくわからんけど，
gcc-c++-2.95.3
では通用したいんちきな書きかたが
gcc-c++-3.3.2
では通用しなくなった
(別の解釈されてしまう)，
というあたりだろうか．
-
gcc-3.3.x
にすると計算結果が違ってくる，
という検索結果はちらほらと見つかる．
しかし，
何が原因なのかよくわからん．
-
どうも leader shoot の開空度がいきなりゼロになる，
というのが原因らしい，
とわかってきた．
ということは，
原因は STL 使いまくりの三次元局所開空度ライブラリの中にアリ，
というような．
いやはや．
難物だ．
-
まぁ，
なンかすぐに計算せんといかんときは
gcc-c++-2.95.3
でやればいいだけなンだけどね．
ということで，
あわててばぐ取りする必要もナシ．
-
時刻はすでに 2030 すぎか．
呪われ言語を夜おそくまででばっぐしてると，
アタマが壊れて生活周期がヘンになっていくので，
本日はここで撤退．
2110 研究室発．
2120 帰宅．
晩飯．
-
気になるので一点だけちぇっく．
局所的明るさを強制的に最大値にしてやると，
それっぽく電気盆栽してんだよなぁ．
ということで，
三次元座標系クラスとかには問題なくて，
たぶん局所開空度計算ライブラリの挙動不審が原因なんだろうなぁ．
[電気盆栽樹齢 24]
-
[今日の運動]
-
エアロバイク 55 分間．
-
腹筋運動 30 ×
3 回．
腕立ふせ 5 ×
3 回．
-
[今日の食卓]
- 朝 (0840):
米麦 0.8 合．
ネギ卵炒飯．
海藻スープ．
- 昼 (1450):
米麦 0.8 合．
ネギ卵炒飯．
海藻スープ．
- 晩 (2240):
パン．
ジャガイモ・タマネギ・ニンジン・ブナシメジ・豆腐のカレー．
2005 年 02 月 07 日 (月)
-
0800 起床．
朝飯．
コーヒー．
0920 自宅発．
曇．
0930 研究室着．
-
三次元開空度ライブラリまわりの調査．
テストコードの再整備 & 試験．
よくわからん．
いったん中断．
時刻は 1115．
[局所開空度計算]
こーゆー単純な試験なんかは問題なくパスするんだがなぁ
……
やはりライブラリ側からではなく，
樹木シミュレイター方面からバグを探索すべきなんだろうか．
しかし，
火急なアレでもないんで，
あとまわしにしますか，
というような．
-
週末にいただいた「おでん」
に関するぎょーむ日誌上の表現に関して，
大学院生からお叱りのメイルいただく．
汗顔のいたりというか
もはや平身低頭にて謝罪状態であります．
これでは、まるで、おいしくないおでんを、くぼさんが残務処理の如く、無理して食べているようではありませんか・・。
ご存知時のとおり、前日から鶏ガラでダシをとり、その上、昆布とカツオのダシもとって（科学調味料は使っていません）、
味付けにも一苦労して作ったんです。みなさんの評価もなかなかだったと思っているのですが。
ちょっとその表現の仕方は、どうかと思いますわ。
えー，
たいへんおいしくいただき，
しかも滋味あふるるものでした．
ゆきとどかぬ表現でいつも皆さんに御迷惑おかけしております．
反省．
それから最後に食べたとき，
具がなかったのは別にとりわけられていたからだそーで
……
うかつにも気づきませんでした．
なンかへんだとは思っていたんですけどねえ
……
-
やはり，
お茶部屋にある食べものには，
食前食後にもつねに細心の注意が必要とされるよな
……
いろいろな意味で．
-
さて，
いつまでも放置しとくわけにはいかんのが，
甲山さんのインドネシアみやげのデータ解析で，
と
……
さっさと計算プログラムとか作ってしまって終らせなくては．
-
あまり進捗せぬまま昼飯．
-
北大生協書籍部にはいつまでたっても到達しなかったので
(たぶん書籍輸送用の馬橇が途中で雪に埋もれてしまったのだろう)
Amazon に発注した統計学本，
「モデル選択
--- 予測・検定・推定の交差点 ---」
(統計科学のフロンティア
シリーズ, 岩波書店)
とどく．
考えてみると，
このしりーづの赤表紙・緑表紙はほとんど全部買ってしまったなぁ．
まあ，
ややまにあっくな内容といってよいので
(ただし説明はとても丁寧である ---
あるていど統計学本に慣れてるヒトたちにとっては)，
そうでもないヒトにはおススめではないかもしれない．
-
3 月の 自由集会，
あまりまにあっくなモノにならなければ良いのだが
……
しかし，
参加者からの基本的な事項の質問が，
しばしばまにあっくな問答になりがちなのである．
前回の自由集会の
質疑応答，
とか
……
これってまさに，
上記「モデル選択」本で下平さんが解説しておられるところだ．
-
うわっ，
美唄方面からきた PDF ファイルひらいたら，
X ごとふりーづしてあせった
……
-
うーむ，
仕事すすまん．
雑事はぼちぼち．
-
仕事ススまぬ理由のひとつは，
アタマの中がまた空間相関推定問題に支配されてるからである．
現時点の私が理解してるこの問題の解決策は 3 つあり
(いづれも Markov Chain Monte Carlo 法を使うのだが)，
-
Bayes 推定
(Maximum a posteriori (MAP) 推定?)
-
Maximum pseudo likelihood estimator 補正による推定
-
EM algorithm を使った推定
(これホントに最尤法になってんのか?)
1. がふつーな方法，
2. が伝統あるけどあまり見かけない方法，
3. はさらに見かけないやつなんだけど，
私のアタマの中で明滅してんのはこれだったりする．
1. はめんどくさそう，
2. はそれに加えて理解が難しい，
というような気がしてですね．
-
Bayes 推定が主流なのはわかるけど，
それだけではイヤだなぁ，
ということで．
-
で，
EM アルゴリズムでこのあたりの問題が解けるとして，
だ．
その中でもいろいろな variant がありうる．
たとえば，
autologistic なモデルのパラメーター推定するとして，
-
Mixed effects
だけを多変量正規分布に拡張
(おもしろくなさそう?)
-
線形予測子全体が
多変量正規分布の乱数
(おもしろそうだけどひどく異端なかんぢがする)
-
二値な隠れ変数が個体数以上に存在し，
これらが推定すべき autologistic model によって生成される
Gibbs 分布にしたがう
(これまたアヤしい --- けど仮定は少なくてすむ)
どれも我流 (は無型?)
なやつらばかりで，
そもそも推定すらマトモにできるかどうかも
計算してみないことにはよくわからんのだけど
……
いちおう，
どれも推定 & 乱数生成のあいだに，
確率論的モデルとしての矛盾はないようにはみえる．
私としては，
3. がヘンてこではあるけれど，
これが一番わかりやすいんでは，
という気がしてですね．
-
東北大の占部さんが甲山さんたちとの相談のためこちらに来ておられる．
お会いするのは初めてで，
ごあいさつする．
占部さんたちが申請しておられる
(そしてなぜかしら私もほんの少しだけまきこまれている)
研究計画の構想について教えていただく．
あれこれ反応がトロい陸上植物群集と異なり，
水の中では何ごとにつけ変動が速いらしい，
ということのよーで．
-
A802 のカラープリンター CANON LBP-2810 の色トナー不調問題，
大学院生と解決に取り組む．
なんかカートリッジ (再生されたもの)
からぼろぼろとこぼれて印刷物に付着してしまうんだよね．
よくわからんけど，
トナーのこぼれかたとカートリッジを観察するに，
シールテイプぬきとりぐちから粉がもれてるよーな気がする．
ガムテイプを切ってそれで封をしてみる
(三色とも)．
はてさて，
どうなることやら．
-
自分のためにも他人のためにもさっさと終わらせたほうがよい，
インドネシアデータ解析のプログラミングの続き．
最尤推定関数をぢりぢりと書く．
しかし，
全体の構想というのがいまいち見えてなくて，
ですね．
-
まだ途中だけど撤退．
1850 研究室発．
1900 帰宅．
運動．
晩飯．
-
[今日の運動]
-
エアロバイク 45 分間．
壊れペダルを修繕したので運動効率が高くなった．
-
[今日の食卓]
- 朝 (0830):
パン．
ジャガイモ・タマネギ・ニンジン・ブナシメジ・豆腐のカレー．
- 昼 (1230):
研究室お茶部屋．
米麦 0.7 合．
トマト・豆腐のスープ．
- 晩 (2210):
スパゲッティー．
ジャガイモ・タマネギ・ニンジン・ブナシメジ・豆腐のカレー．
2005 年 02 月 08 日 (火)
-
0700 起床．
朝飯．
コーヒー．
また MCMC 問題の検討に没頭してしまう．
0920 自宅発．
曇．
気温高い．
0930 研究室着．
-
インドネシアデータ解析用のプログラム書き
in R．
1150 どーにかこーにか中核部分というか
いちばん面倒そうなあたりできた．
あとは，
少しひねるだけでいいんだろう．
たぶん．
今日は甲山さんきてないんで，
これ以上はススまん．
[いまいちアヤしい]
直径成長の分布．
まあ，
全 159 樹種こみだから，
なのかもしれないけど．
直径のばらつきはガンマ分布を仮定．
平均は exp(線形予測子)．
offset 項として観測間隔を指定
(というのはちょっとイヤなんだけど，
割り算値とか使うよりはマシなんで)．
分散は平均の定数倍
(このデータでは分散は平均より小さい)．
よって，
glm()
とか使えない．
こういう場合は尤度関数を自分で書き，
optim()
など使って最尤推定値を計算する．
2005 年 02 月 09 日 (水)
-
0830 起床．
朝飯．
コーヒー．
0930 自宅発．
曇．
0940 研究室着．
-
昨晩はまりこんでしまった混合モデル推定問題のつづき．
どうも個体 in ブロックという二階層 random effects
はなかなかうまく推定できんようなので，
とりあえずブロックの効果だけのモデルに単純化してみる．
-
random effects をひとつにすれば，
最尤法を使った混合モデルの推定ができる．
乱数を使ってブロック差・個体差のあるデータを生成する．
ブロック (plot) のランダム効果は
rnorm(1, mean = 0, sd = 0.5)
で，
個体のそれは
rnorm(1, mean = 0, sd = 0.1)
という設定にしている．
いつものことながら，
この同じデータを使っていても
glmm() (repeated package)
の推定結果と
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.586 0.258 2.27 0.0231
toxicity.L 3.646 0.653 5.58 2.4e-08
toxicity.Q -0.840 0.612 -1.37 0.1701
toxicity.C -1.482 0.635 -2.33 0.0196
toxicity^4 1.176 0.642 1.83 0.0669
toxicity^5 0.306 0.635 0.48 0.6301
sd 0.764 0.266 2.87 0.0041
と glmmML() (glmmML package)
の推定結果
coef se(coef) z Pr(>|z|)
(Intercept) 0.580 0.291 1.992 4.6e-02
toxicity.L 3.614 0.822 4.396 1.1e-05
toxicity.Q -0.834 0.679 -1.230 2.2e-01
toxicity.C -1.471 0.710 -2.071 3.8e-02
toxicity^4 1.168 0.709 1.648 9.9e-02
toxicity^5 0.305 0.695 0.439 6.6e-01
Standard deviation in mixing distribution: 0.726
はびみょーに異なる
(AIC はだいたい同じ)．
glmmPQL() (MASS package)
だとこんなかんぢで．
Random effects:
Formula: ~1 | plot
(Intercept) Residual
StdDev: 1.0456 0.81762
Fixed effects: dead ~ 1 + toxicity
Value Std.Error DF t-value p-value
(Intercept) 0.5811 0.27448 84 2.1173 0.0372
toxicity.L 3.5539 0.66403 36 5.3520 0.0000
toxicity.Q -0.7927 0.66284 36 -1.1959 0.2395
toxicity.C -1.4142 0.67645 36 -2.0907 0.0437
toxicity^4 1.1011 0.68171 36 1.6152 0.1150
toxicity^5 0.2699 0.67646 36 0.3989 0.6923
-
同じデータをもちいて，
こんどはブロック (plot) 差だけでなく個体差 (id) も考慮して
glmmPQL()
で推定させると
(推定のさせかたは昨日のぎょーむ日誌
参照)，
Random effects:
Formula: ~1 | plot
(Intercept)
StdDev: 5.9131e-06
Formula: ~1 | id %in% plot
(Intercept) Residual
StdDev: 6.7657 4.4077e-06
Fixed effects: dead ~ 1 + toxicity
Value Std.Error DF t-value p-value
(Intercept) 1.5665 0.61762 84 2.5364 0.0131
toxicity.L 12.9624 1.51286 36 8.5681 0.0000
toxicity.Q -2.0533 1.51286 36 -1.3572 0.1832
toxicity.C -4.1781 1.51286 36 -2.7618 0.0090
toxicity^4 3.0496 1.51286 36 2.0158 0.0513
toxicity^5 -0.8537 1.51286 36 -0.5643 0.5761
このように，
むちゃくちゃになる
(とくに random effects)．
-
これはいくら何でもひどい
……
ということで
glmmPQL()
の使いかたを間違ってる可能性を検討してみる
……
うーむ，
ヘンなことはやってないハズ，
なんだが．
しかし，
ブロック内個体数を多くする設定のもとでも，
事態は改善されんのだよねえ．
(後記:
これは glmmPQL()
の使いかたが間違ってるとゆーより，
私のモデルのたてかたが良くなかった，
とわかった
……
これでは推定できなくて当然
……
詳しくはこの週末のぎょーむ日誌)
-
さらに，
だ
……
上のいずれの推定においても fixed effects
の推定値もよくわからないものになっている．
これは toxicity なる処理を順序因子
(
ordered)
として表現してる影響のよーで
……
たとえば (Intercept)
をここに書いてあるとーり，
と理解してはいかん．
うーむ，
選点な直交多項式 (orthogonal polynomial)
め．
-
ナゾに満ちた選点直交多項式系，
John Fox ``An R and S-PLUS Companion to Applied Regression''
によると，
これが使える必要条件は
-
水準間が「等間隔」
-
各水準における標本数が同じ
だそーだ．
1. は知っていただけど，
2. は知らなかった．
嗚呼，
苫小牧某修論とかどうなることやら
……
ともあれ，
R
における順序因子の取りあつかいは，
いっそう慎重でなきゃならん，
とわかった．
このあたり，
単純ではないよーで，
慶應大の統計学講義ノート PDF ファイルには
「あえて非順序因子として」
うんぬんと書かれていたり，
あるいは John Fox は ``Helmert contrasts may also be of interest for
an ordered factor.'' などと注記している
---
R では Helmert 対比は非順序因子 (factor() な連中)
の default の coding として使われてるアレ．
-
とはいえ，
いま乱数例を作って調べてる問題は，
上記の 1. も 2. も完ペキに満たしてるんだけどねえ
……
だれかさっさと
biometry
に模範解答を投稿してくれ．
-
たしかに Helmert 対比のほうがわかりやすい
……
しかし，
矛盾する結果が出やすいのもこちらだろう．
-
三浦さんの ThinkPad に Norton AntiVirus
インストールする作業．
二回も再起動を要求された．
あいかわらず，
ゐんどーづはヘボい．
-
研究費の残り 2800 円，
か．
いつのまにか消失してますなぁ．
昼飯．
-
昨日の私の軌道おんちぶりを見かねて，
プラネテス
でも鑑賞して勉強しなさい，
という助言メイルいただく
……
とある女性大学院生から．
これが理系わーるどだ．
-
インドネシアなデータ解析の下請けぎょーむ．
発注者側，
なんかだんだん意味不明になってきた．
とりあえず，
計算プログラムまわりの作業は終了．
あとは引き渡して終りにしますか．
とりあえず，
1740 当方の作業は終了．
ここからは撤退．
-
……
いや，
待て．
この計算に関する文書かき作業が残されていた，
か
……
作業終了 1810．
オワった．
-
1850 研究室発．
1900 帰宅．
運動．
晩飯．
明日の輪読会の予習．
-
お．
粕谷さんが biometry に投稿．
しかし現時点では指摘 3 点のみ，
か．
-
GLM こそが不等分散な状況をあつかう
-
「不等分散ならのんぱら」信仰はまちがい
(cf. 先月末のぎょーむ日誌)
-
生存時間のモデル化したいんだったら負の指数分布とか?
……
負の指数分布ってあったのか
……
どうも調べてみると，
ふつーの指数分布をそう呼んでる場合もあれば，
生存時間解析で使ってるあの分布
(私が勝手に「二重な指数分布」と読んでるやつ)
のことのよーで．
しかし，
おそらくこの問題は時間分布ではあるまい．
(後記:
またしても，
まちがい
……
時間分布のハナシである)
-
[今日の運動]
-
[今日の食卓]
- 朝 (0850):
米麦 0.5 合．
卵焼．
- 昼 (1350):
研究室お茶部屋．
米麦 0.6 合．
ジャガイモ・ニンジン・タマネギ・大豆のシチュー．
- 晩 (2130):
米麦 0.7 合．
ジャガイモ・ニンジン・タマネギ・大豆のシチュー．
2005 年 02 月 10 日 (木)
-
0720 起床．
朝飯．
コーヒー．
0830 自宅発．
晴．
0840 研究室着．
輪読会の予習．
-
1000
Life History Evolution in Plants
輪読会．
本日の担当は M1 実吉さんで Chapter 7.
Evolution of plant dispersal．
大陸-島モデルの説明にでてくる mainland の訳語が
「内地」
であるところが北海道的だ．
-
1255 終了．
昼飯．
雑書類処理少々．
いかん，
輪読会ばてぎみだ．
-
大学院生のゐんどーづ機でプリンター設定．
EPSON のヘンてこプロトコルを TCP/IP 経由でやりとりする方法がわからん．
めんどうになったので，
プリンター横のゐんどーづ機をケイブルで直結．
「共有」とやらを設定．
他のホストからは
この「共有」されたプリンターを通して印刷するようにしてみる．
うーむ
……
-
いきなりどか雪．
窓の外の視界が閉ざされる．
-
またくだんの EPSON LP-2000C とらぶる
……
印刷できないので現況報告させてみると，
「廃トナーカートリッジ」とやらが満杯状態とのこと．
塵ぱいと粉塵爆発の危険におびえつつ，
廃棄処理．
一件落着．
-
統計学勉強．
EM algorithm まわり．
-
1900 ごろ香川大の小林さんと低温研の加藤さんみえる．
とりあえず晩飯にでる．
外は吹雪．
北 13 西 4 の瀋陽飯店．
低温研の飯村さんと合流して晩飯．
-
また吹雪のなか地環研までもどって，
研究室お茶部屋にて母子里とアカマツ相談．
母子里に関しては Monte Carlo Markov Chain (MCMC) 法でもって，
葉配置の Gibbs 分布を生成する，
というアイデアを説明する．
計算やるまえからすでに「手ぬき」のやりかたあれこれ検討．
アカマツに関しては分割整理案の確認，
というような．
-
雑談してたらすっかり遅くなってしまった
……
2640 研究室発．
降雪，
すこし弱まる．
2700 帰宅．
-
[今日の運動]
-
[今日の食卓]
- 朝 (0740):
米麦 0.7 合．
ジャガイモ・ニンジン・タマネギ・大豆のシチュー．
- 昼 (1310):
研究室お茶部屋．
米麦 0.7 合．
ジャガイモ・ニンジン・タマネギ・大豆のシチュー．
- 晩 (1950):
北 13 西 4 の瀋陽飯店．
ニラレバ定食 680 円．