2012年12月17日

たいてい土曜日は出勤して勤務日にできない集中して長時間できる仕事に取り掛かることにしています。というのも、周りに人はいないし、電話も上司からの依頼もないから。ただ、昼過ぎから部下が2人ほど出てきて仕事をしていたので、まるっきりの1人というわけではなかった。

今日は、レズリー行列のプログラムにとりかかる。共同研究者から送られてきたプログラムは最適化はされているのだが、どうも微妙に結果が思うように得られない。あれこれと、コードを抽出して走らせていると、根本のレズリー行列式自体が間違っていたという落ち。さすがに間違っていないだろうという思い込みがまずかった。ここを修正して、いくつか気になるところを触っていたら、あっというまに数時間ほど経っていました。

レズリー行列は生態学などでは結構使うと思うのですが、Rではあまり紹介されているブログなどがない。Rjpwikiにも多少あるが、ちょっと長ったらしい。こんど私が使っているコードも載せてみようかと思っているが、野生生物の個体群動態に関するコードも少しずつ書いてみよう。
まずは超簡単なところから。とある野生生物の個体群があり、それが増加率αによって増加しているパターン。

rm(list=ls())
period <- 100
rate <- 1.1
initial.size <- 20
pop.size <- numeric(period)
for (i in 1:period){
pop.size[1] <- initial.size
pop.size[i+1] <- pop.size[i]*rate
}
plot(pop.size, xlab="period", type="l")

これで下のような図になるはず。

これは増加率1.1/世代時間で、20世代時間ほど動かした例。普通、こんな単調増加は示さない。増加率は環境要因などによって変わり得る。したがって、増加率そのものをランダムに0.3〜2まで発生すると仮定する。
rm(list=ls())
set.seed(10)
period <- 100
initial.size <- 20
max.rate <- 2
rate <- runif(period,min=0.3,max=max.rate)
generation <- seq(from=1, to=period)
Npop <- matrix(0,period,1)
Npop <- initial.size
for (i in 2:period){
Npop[i] <- rate[i-1]*Npop[i-1]
}
plot(generation, Npop, xlab="generation",type="l")
これで下のようになるはず。

set.seed()の括弧内の数字をかえることで、色々なパターンが試せるはず。

今日はこんなところで。特に難しいコードではありません。

朝：トースト、昼：近所のパン屋さんで仕入れたサンドイッチなど、夜：カレーライス

勉強：なぜかみんなテレビをみていて、5分ほど4年生の子供に算数をさせたら半なきだった。どうして、5分くらいの勉強で泣く！？