メモ。統計ソフトRを使って、個体群の絶滅確率を計算する関数を考え中。
n回繰り返しを行って、そのうち絶滅した回数を数えればよい。
f <- function (n) {
count <- 0
for (i in 1:n) { # 1~nの値をとりながら{}内を繰り返す
if (絶滅する条件) count <- count +1
}
return (count/n)
}
肝心の「絶滅する条件」をどう書けばよいのか、まだ分かりません。
（追記：2015.02.02）
行列の累乗の書き方。固有ベクトルの関数 eigen()$vecotrsと、固有値の対角行列diag(eigen()$values)、固有ベクトルの逆行列 solve(eigen()$vectors)を用いて下のように書けます。
A^2 <- eigen(A)$vecotrs %*% diag(eigen(A)$values**2) %*% solve (eigen(A)$vectors)
てことで初期の齢分布をwベクトル、個体群の時間変化を示す射影行列をA（乱数による変動の項も組み込んでおく）とすると、
f <- function (n) {
count <- 0
A^n <- eigen(A)$vecotrs %*% diag(eigen(A)$values**n) %*% solve (eigen(A)$vectors)
for (i in 1:n) {
if (A^n %*% w == 0) count <- count +1
}
return (count/n)
}
（追記：2015.02.14）
行列の要素の置換。
> A<-matrix(5, nrow=3, ncol=3)
> A
[,1] [,2] [,3]
[1,] 5 5 5
[2,] 5 5 5
[3,] 5 5 5
> A[2,2]<-c(6) # 2行2列目だけを置換。
> A
[,1] [,2] [,3]
[1,] 5 5 5
[2,] 5 6 5
[3,] 5 5 5
射影行列Aにファットヘッドミノーの生命表を使う（Miler and Ankley, 2004）。行列Aの要素[1,1]だけを正規分布に従う乱数に変更してみる。平均値は0.75のまま。1000年分計算。
if関数からifelse関数に変更。if関数のままだとベクトルが定義できなくてエラーが出る。
f <- function (n) {
yr<-1000
count <- 0
A <- matrix(c(0.75,0.39,0, 1.5,0,0.39, 3,0,0), nrow=3, ncol=3)
w <- eigen(A)$vectors [,1]
Ayr <- eigen(A)$vectors %*% diag(eigen(A)$values**yr) %*% solve (eigen(A)$vectors)
for (i in 1:n) {
A[1,1] <- c(rnorm(1, mean=0.1, sd= 0.02))
A[1,1] <- c(rnorm(1, mean=0.1, sd= 0.02)) # 非負にするためsd設定
ifelse (Ayr %*% w == 0, count <- count +1, count <- count)
}
return (count/n)
}
（追記：2015.02.15）
1000年分計算するならば、個体群の大きさの時系列変動をグラフにしてみたい。
A <- matrix(c(rnorm (1, mean=0.75, sd= 0.2),0.39,0, 1.5,0,0.30, 3,0,0), nrow=3, ncol=3)
w <- eigen(A)$vectors [,1]
number <- function (x)
Re(colSums(eigen(A)$vectors %*% diag(eigen(A)$values**x) %*% solve (eigen(A)$vectors) %*% w))
plot (number, 0, 1000)
number()ではちゃんと計算できるのに、なぜかplot(number,,)にすると下のエラーが出る。
以下にエラー eigen(A)$vectors %*% diag(eigen(A)$values^x) :
適切な引数ではありません
追加情報: 警告メッセージ:
In eigen(A)$values^x :
長いオブジェクトの長さが短いオブジェクトの長さの倍数
（追記：2015.02.17）
良い例を見つけた。神戸大学の中澤港教授という先生のサイト。マラリア感染率のシミュレーションモデルをRで作成しています。これを参考にして個体群サイズの時間変化シミュレーションをおこなってみます。
（追記：2015.02.18）
上の中澤先生のコードを参考にして改良してみました。個体群サイズnumberをシミュレーション年数分のベクトルとして定義しないとダメなんですね。
number <- function (yr) {
N <- rep (0,yr)
A <- matrix(c(0.75,0.39,0, 1.5,0,0.30, 3,0,0), nrow=3, ncol=3)
w <- floor(5*Re(eigen(A)$vectors [,1])) # 初期年齢構成
N[1] <- sum(w) # 1年目の個体数
for (i in 2: yr) {
a <- rnorm (1,mean=0.75, sd= 0.25)
A <- matrix(c(a,0.39,0, 1.5,0,0.30, 3,0,0), nrow=3, ncol=3)
w <- A%*%w
N[i] <- sum(w)
}
N <- floor(N)
}
sim <- function (n, yr) {
No <- number (yr)
plot (1:yr, No, type="l", xlab="year", ylab="Number of individuals")
for (j in 2:n) {
No <- number (yr)
points (1:yr, No, type="l")
}
}
シミュレーション結果の例。10年分を50回やってみました。
> sim(50,10)
ただ今気づきましたが、これまで書いてきたプログラムだと個体数が負になっても計算し続けてしまいます。例えば 「 a <- rnorm (1,mean=0.75, sd= 1.0)」のように偏差を大きくして計算した下の結果。IF関数を導入して修正しなければ。
> sim(5,30)
（追記：2015.02.19）
上に書いた反省点を修正しました。なんかまどろっこしいコードになったけど、まあいいや。
number <- function (yr) {
N <- rep (0,yr)
A <- matrix(c(0.75,0.39,0, 1.5,0,0.30, 3,0,0), nrow=3, ncol=3)
w <- floor(10*Re(eigen(A)$vectors [,1])) # 初期年齢構成
N[1] <- sum(w) # 1年目の個体数
for (i in 2: yr) {
a <- rnorm (1,mean=0.75, sd= 1.0)
A <- matrix(c(a,0.39,0, 1.5,0,0.30, 3,0,0), nrow=3, ncol=3)
w <- floor(A%*%w)
ifelse (w <=0, 0 , w) # wの0以下の要素を0にする
N[i] <- floor(sum(w))
if (N[i] <=0) {
N[i:yr] <- 0
w[1:nrow(w)] <-0
} else {
N[i]
}
}
N
}
sim <- function (n, yr) {
No <- number (yr)
plot (1:yr, No, type="l", xlab="year", ylab="Number of individuals", ylim = c(-1, max(No)))
for (j in 2:n) {
No <- number (yr)
points (1:yr, No, type="l", ylim = c(0, max(No)))
}
}
（追記：2015.02.20）
上の修正点を、絶滅確率のカウント関数にも適用しました。2/19追記分のnumber関数をそのまま使う。
prob <- function (n, yr) {
count <-0
for (j in 1:n) {
finalnumber <- number(yr)
ifelse ( finalnumber == 0, count <- count+1, count <- count)
}
return (count/n)
}
plot（ここではsimという関数）と上記のprob関数の結果が合致しているかどうか検証してみました。
> prob(10000,5)
[1] 0.3201
まずprob関数の結果。5年間での絶滅確率は10000回シミュレーションした結果、およそ0.3でした。プロットでも同じくらいの確率になっています。シミュレーション回数を多くすると線が重なってしまうので、10回のシミュレーションを何回かおこなってざっくり確認してみました。
> sim(10,5)
> sim(10,5)
