1 R 事はじめ

1.1 Pella-Tomlinson 余剰生産モデル (確率変動なし)

1.1.1 モデルの定義

\[ P_{t+1} = P_t + r P_t \left\{1-\left(\frac{P_t}{K}\right)^z \right\} - C_t \]

1.1.2 初歩的に漸化式を利用して個体群動態を計算

r <- 0.2
K <- 10000 
z <- 1
TT <- 50

Catch <- rep(500, TT)

TT1 <- TT-1
P <- numeric(TT)
P[1] <- K
for(t in 1:TT1)
{
 tmp <- P[t] + r*P[t]*(1-(P[t]/K)^z) - Catch[t] 
 P[t+1] <- max(tmp, 0.001)
}

print(P, digits=0)
##  [1]10000 9500 9095 8760 8477 8235 8026 7843 7681 7537 7409 7293 7187 7092
## [15] 7004 6924 6850 6781 6718 6659 6604 6552 6504 6459 6416 6376 6338 6303
## [29] 6269 6236 6206 6177 6149 6123 6097 6073 6050 6028 6007 5987 5967 5949
## [43] 5931 5913 5897 5881 5865 5850 5836 5822

1.1.3 計算結果の図示

plot(1:TT, P)

plot(1:TT, P, type="l", lwd=2, col="blue", xlab="Year", ylab="Population size", ylim=c(0,K))

1.1.4 もう少し進んだ方法：個体群動態計算用関数PDM.PTを作成

#Pella-Tomlinson Production model (without stochastic process error)
PDM.PT<-function(r, K, z, P1, TT, Catch)
{
 TT1 <- TT-1
 P <- numeric(TT)
 P[1] <- P1
 for(t in 1:TT1)
 {
  tmp <- P[t] + r*P[t]*(1-(P[t]/K)^z) - Catch[t] 
  P[t+1] <- max(tmp, 0.001)
 }
 return(P)
}

1.1.5 新しく作った個体群動態計算用関数PDM.PTを実行

Res <- PDM.PT(r=0.2, K=10^4, z=1, P1=10^4, TT=50, Catch=rep(500,TT))
print(Res, digits=0)
##  [1]10000 9500 9095 8760 8477 8235 8026 7843 7681 7537 7409 7293 7187 7092
## [15] 7004 6924 6850 6781 6718 6659 6604 6552 6504 6459 6416 6376 6338 6303
## [29] 6269 6236 6206 6177 6149 6123 6097 6073 6050 6028 6007 5987 5967 5949
## [43] 5931 5913 5897 5881 5865 5850 5836 5822
plot(1:TT, Res, type="l", lwd=2, col="blue", xlab="Year", ylab="Population size", ylim=c(0,K))

1.2 Pella-Tomlinson 余剰生産モデル (確率変動あり)

1.2.1 モデルの定義

\[ P_{t+1} = \left[ P_t + r P_t \left\{1-\left(\frac{P_t}{K}\right)^z \right\} - C_t \right]*\exp(\varepsilon_t), \quad \varepsilon_t \sim N(0, \sigma^2) \]

1.2.2 さっき作った個体群動態計算用関数PDM.PTを書き換え（上書き更新）

#Pella-Tomlinson Production model (with stochastic process error)
PDM.PT<-function(r, K, z, P1, TT, Catch, sigma=0)
{
 TT1 <- TT-1
 P <- numeric(TT)
 Epsilon <- rnorm(n=TT, mean=0, sd=sigma)
 P[1] <- P1
 for(t in 1:TT1)
 {
  tmp <- P[t] + r*P[t]*(1-(P[t]/K)^z) - Catch[t] 
  tmp <- tmp*exp(Epsilon[t])
  P[t+1] <- max(tmp, 0.001)
 }
 return(P)
}

1.2.3 再定義したPDM.PTを実行

K <- 10^4
TT <- 50

Res <- PDM.PT(r=0.2, K=K, z=1, P1=K, TT=50, Catch=rep(500,TT), sigma=0.05)
print(Res, digits=0)
##  [1]10000 8906 9231 9374 9148 8889 7749 6752 6759 6727 5945 5965 6844 6678
## [15] 6922 6562 6463 6270 6149 6339 6452 6924 6582 6572 5906 6569 6489 7041
## [29] 7071 6705 6207 5996 6002 6302 5994 6052 5577 5752 5333 5262 5192 5075
## [43] 4927 5088 5240 6205 6160 6046 6427 5703
plot(1:TT, Res, type="l", col="blue", xlab="Year", ylab="Population size", ylim=c(0,K))

Res <- PDM.PT(r=0.2, K=K, z=1, P1=K, TT=50, Catch=rep(500,TT), sigma=0.05)
print(Res, digits=0)
##  [1]10000 8732 8541 9107 9041 7880 7590 7943 7515 7139 6946 6582 5713 5536
## [15] 5442 5499 5572 5657 5222 5407 5538 5511 6135 6693 6306 6475 6707 6481
## [29] 5872 5755 5819 6021 6142 6629 6777 6543 6331 6540 6546 6215 5757 5486
## [43] 5410 5717 5606 6206 6643 6202 5772 5436
points(1:TT, Res, type="l", col="blue")

1.2.4 PDM.PTを繰り返し実行して結果を図示

Nsim <- 100
Pmat <- array(NA, c(Nsim, TT))
plot(0, type="n", xlab="Year", ylab="Population size", xlim=c(0,TT),  ylim=c(0,K))

for(i in 1:Nsim)
{
 Pmat[i,] <- PDM.PT(0.2, 10^4, 1, 10^4, 50, rep(500,TT), 0.05)
 points(1:TT, Pmat[i,], type="l", col="lightblue")
}
P.med <- apply(Pmat, 2, median); 
P.L5per <- apply(Pmat, 2, quantile, 0.05)
P.U5per <- apply(Pmat, 2, quantile, 0.95)
points(1:TT, P.med, type="l", lwd=2, col="red")
points(1:TT, P.L5per, type="l", lty=2, lwd=2, col="red")
points(1:TT, P.U5per, type="l", lty=2, lwd=2, col="red")

polygon(
  c(1:TT, TT:1), c(P.L5per, rev(P.U5per)), 
col = "#00000020", border = NA)

1.2.5 上記を関数化

PDM.PT.sim <- function(r, K, z, P1, TT, Catch, sigma=0, Nsim)
{
 Pmat <- array(NA, c(Nsim, TT))
 plot(0, type="n", xlab="Year", ylab="Population size", 
      xlim=c(0,TT), ylim=c(0,K), xaxs="i", yaxs="i")

 for(i in 1:Nsim)
 {
  Pmat[i,] <- PDM.PT(0.2, 10^4, 1, 10^4, 50, rep(500,TT), 0.05)
  points(1:TT, Pmat[i,], type="l", col="lightblue")
 }
 P.med <- apply(Pmat, 2, median); 
 P.L5per <- apply(Pmat, 2, quantile, 0.05)
 P.U5per <- apply(Pmat, 2, quantile, 0.95)
 points(1:TT, P.med, type="l", lwd=2, col="red")
 points(1:TT, P.L5per, type="l", lty=2, lwd=2, col="red")
 points(1:TT, P.U5per, type="l", lty=2, lwd=2, col="red")

 polygon( c(1:TT, TT:1), c(P.L5per, rev(P.U5per)), col = "#00000020", border = NA)
 
 rect(xleft=-1, ybottom=0, xright=TT, ytop=0.3*K, lwd=0, col=rgb(1,0,0,alpha=0.2))
 rect(xleft=-1, ybottom=0.3*K, xright=TT, ytop=0.5*K, lwd=0, col=rgb(1,1,0,alpha=0.2))
 rect(xleft=-1, ybottom=0.5*K, xright=TT, ytop=K, lwd=0, col=rgb(0,1,0,alpha=0.2))
 
 return(Pmat)
}

1.2.6 実行！

Res <- PDM.PT.sim(0.2, 10^4, 1, 10^4, 50, rep(500,TT), 0.05, Nsim=1000)