備忘録 a record of inner life

やったことや考えたこと・本・論文・音楽の備忘録。 特に環境科学・生態毒性に関して。

Rを用いた絶滅確率のモンテカルロシミュレーション:その2

Rで遊ぶのがすっかり楽しくなってきました。

 

その1に続いて、ファットヘッドミノーの生命表を使います(Miller and Ankley, 2004)。今度は初期個体数による絶滅確率の違いを検討してみます。その1ではw(9,2,0)の初期個体数11匹で計算していましたが、その年齢構成比率は変えずに数だけを変更するとどうなるでしょう。

確率分布は1~3齢目の産仔数に組み込みます(下の太字の部分)。初期個体数はxで調整できる関数として設定。

number  <- function (x, 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(x*Re(eigen(A)$vectors [,1]))   

                  N[1] <- sum(w)            # 1年目の個体数

     for (i in 2: yr) {

                        a <- rnorm (1,mean=0.75, sd= 0.15)                                                  

                        b <- rnorm (1,mean=1.5, sd= 0.3)

                        c <- rnorm (1,mean=3, sd= 0.6)

                        A <- matrix(c(a,0.39,0, b,0,0.30, c,0,0), nrow=3, ncol=3)

                        w <- floor(A%*%w) 

                        ifelse (w <=0, 0 , w)  

                        N[i] <- floor(sum(w))

                        if (N[i] <=0) {

                             N[i:yr] <- 0

                             w[1:nrow(w)] <-0

                           } else {

                            N[i]

                           }

                       }

        N

}

sim  <- function (x, n, yr) {

   No <- number (x, 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 (x, yr)

                       points (1:yr, No, type="l", ylim = c(0, max(No)))

          }

}

prob <- function (x, n, yr) {

   count <-0

            for (j in 1:n) {

                       finalnumber <- number(x, yr)

                       ifelse ( finalnumber == 0, count <- count+1, count <- count)

          }

 return (count/n)

}

 

結果。25年間の絶滅確率を計算してみました。

今回の計算条件では、初期個体数が5匹(下のprob(5,1000, 25))だと25年の間にファットヘッドミノーが絶滅する確率は約30%でした。初期個体数が増えるにつれ、絶滅確率も減少し、1年目に10~11匹もいれば25年後の絶滅確率はほぼ0になる、と計算できました。

> prob(5,1000,25)

[1] 0.336

> prob(6,1000,25)

[1] 0.233

> prob(7,1000,25)

[1] 0.053

> prob(8,1000,25)

[1] 0.006

> prob(9,1000,25)

[1] 0.001

> prob(10,1000,25)

[1] 0

> sim(10,1000,5)

f:id:Kyoshiro1225:20150223084930j:plain

 

 

(追記 2015.02.26)

上のコードでは

                      w <- floor(A%*%w) 

                        ifelse (w <=0, 0 , w)  

                        N[i] <- floor(sum(w))

 と書いてましたが、これだと過小評価かも。2回もfloor関数での切り捨てをおこなっているので。下のように修正すべきですね。

                      w <- A%*%w

                        ifelse (w <=0, 0 , w)  

                        N[i] <- floor(sum(w))

 このコードだと、上記の計算条件でも初期個体数が3匹さえいれば絶滅確率は0%になりました。