2014年2月1日

晴れ(すこし暖かい気がする)

朝、やっぱり土日が起きられない。それでも6時50分には目が覚めて、朝食をとり身支度。来週木曜日から家に戻って休息するため、とりあえずこの週末も仕事をする。8時前に職場について、今月のWSのプレゼンテーションのたたき台を作成し、共同研究者に送る。つぎに、先月のアメリカでの会議報告書を作成し、一緒に出張したメンバーに送りました。この時点で11時だったので、いったんやめて、お弁当を買いに。暖かく気持ちよかったので、あるいて10分ほどのスーパーに行きました(いつもは、自転車)。

昼食後、Leslie行列を用いたシミュレーションのプログラムを書くためにあれこれテストコードを。

先日購入した、

生物数学入門 ?差分方程式・微分方程式の基礎からのアプローチ?

生物数学入門 ?差分方程式・微分方程式の基礎からのアプローチ?

行列式の組み立て方の参考になります。

一回、Leslie行列のプログラムをRで書いたことはあるのですが、その時は行列の中身を1つ1つ書いています。寿命がせいぜい10歳くらいだからそのようなことが出来ましたが、今回の生物は寿命が30歳以上。ということで、すこし行列の中身の書き方を工夫してみました。


  L <- matrix(0,ncol=longevity, nrow=longevity) # initical Leslie matrix with non uncertainty
# step1 Leslie matrix with survivorship
for (j in 1:longevity){
S[j] <- S*unc.survive[j]
}
for (i in 1:longevity){
L[i+1,i] <- exp(-S[i])
}

いったん、Lという行列の中に、0をつめてしまって、Sの生残式をLeslie行列の中に詰め込んでしまいます。unc.surviveはなんでもいいのですが、誤差項です。いまのところ対数正規誤差を使ってます。こんな感じにすると、年齢毎の生産率が誤差項を持って表すことが可能です。次は、加入部分の行列を同じように作って、先ほどの行列と足し算をすることでLeslie行列の完成です。すぐに忘れてしまうので、メモっときました。プラスグループにしたいときには、最終年齢のところに同じ値を詰め込めば良いです。このところは、あとでシナリオとして付け足しておきましょう。

いま見たいのは、生産率と加入量が誤差を伴っているときに、内的自然増加率がどの程度の幅でばらつくのか?をシミュレーションするだけです。今のところのアイデアは、10000年くらい回してみて、初期の1000年を除く9000年で内的自然増加率を拾ってくれば良いのでは?と考えていますが、ちょっとオーバーメモリする可能性もありますね。初期の年齢組成を複数与えておいても、たぶん100年くらいで安定齢に近くなると思っています。
テスト結果

なにかプログラムミスしてそう・・・。100年くらいで増加率は一定になっていくようだけど、気持ち悪い。

ただし、これだと計算期間中の行列内の数値は一緒になるので、それを変えるように工夫してみよう。

とりあえず、2月の上旬の目標をプログラム作成にしよう。

Springer本はY先生にとりあえず原稿を送りました。ちょっと肩の荷がおりました。

 朝:トースト
 昼:焼きそば
 夜:鯖の刺し身

勉強:EXオープン模試(6年生と3年生)。大手新聞は首都圏の中学受験の報道が目立ちますね。我が家はすでに終わって進学先も決まっているので、気持ちが楽です。終わってみれば、受験をすることで得られることもあると思います。もちろん、受験しなくても義務教育で公立中学校で一定の教育を受けることは可能です。