★新サイト完成しました!
３秒後に自動的に移動します
変わらない方は こちらからどうぞ
http://logics-of-blue.com/vpa%EF%BC%88%E5%8C%97%E5%A4%A7%E4%BB%95%E6%A7%98%EF%BC%89/
資源量推定手法の王道、ＶＰＡのプログラムです。
最終齢はプラスグループで漁業が定常な場合を想定しました。プラスグループの計算においてはPopeの近似式は使って いません。
最近年のFは過去3年間のFの平均値だと仮定。ま た、最高齢のFは最高齢−1齢のFと等しいと仮定しています。
#===========================
# VPA プログラム
#===========================
VPA<-function(data,M,F.first){
data<-data
M<-M
F.first<-F.first
#まずは普通に計算
VPA.first<-function(data,M,F.first){
C<-data
#漁獲量
N<-data
#資源尾数
F<-data
#漁獲係数
F.first<-F.first
#最新年・最高齢魚のFの設定
F[nrow(F),ncol(F)]<-F.first
M<-M
#自然死亡係数
for (j in 0:(ncol(C)-2)){
N[nrow(N),ncol(N)-j]<-C[nrow(C),ncol(C)-j]*(F[nrow(F),ncol(F)-j]+M)/F[nrow(F),ncol(F)-j]
# ↑ まずは最高齢魚のNの計算
for (i in
1:nrow(C)){
#
↓ シングルコホート解析
if ((ncol(C)-i-j)<=1)break
N[nrow(N)-i,ncol(N)-i-j]<-N[nrow(N)-i+1,ncol(N)-i+1-j]*exp(M)+C[nrow(C)-i,ncol(C)-i-j]*exp(M/2)
}
for (i in
1:nrow(C)){
# ↓ Fの計算
if ((ncol(C)-i-j)<=1)break
F[nrow(F)-i,ncol(F)-i-j]<-log(N[nrow(N)-i,ncol(N)-i-j]/N[nrow(N)-i+1,ncol(N)-i+1-j])-M
}
if
((ncol(C)-j)<=2)break
F[nrow(F),ncol(F)-j-1]<-F[nrow(F)-1,ncol(F)-j-1]
# ↑ Fの近似
}
for (j in (1:(nrow(C)-1))){
F[nrow(F)-j,ncol(F)]<-mean(c(F[nrow(F)-j,ncol(F)-1],F[nrow(F)-j,ncol(F)-2],F[nrow(F)-j,ncol(F)-3]))
# ↑ 最新年Fを、3年間の平均値と近似
N[nrow(N)-j,ncol(N)]<-
C[nrow(C)-j,ncol(C)]*(F[nrow(F)-j,ncol(F)]+M)/(F[nrow(F)-j,ncol(F)]*(1-exp((-1)*(F[nrow(F)-j,ncol(F)]+M))))
# ↑ 近似したFから、最新年のNの推定
if
((nrow(C)-1-j)==0) break
for (i in
1:(nrow(C)-1-j)){
#
↓ シングルコホート解析
N[nrow(N)-j-i,ncol(N)-i]<-N[nrow(N)-j-i+1,ncol(N)-i+1]*exp(M)+C[nrow(C)-i-j,ncol(C)-i]*exp(M/2)
}
for (i in
1:(nrow(C)-1-j)){
# ↓ Fの計算
F[nrow(F)-i-j,ncol(F)-i]<-log(N[nrow(N)-i-j,ncol(N)-i]/N[nrow(N)-i+1-j,ncol(N)-i+1])-M
}
}
A<-list(推定資源尾数=N,推定漁獲係数=F)
return(A)
}
#最適化のための、残差平方和を計算
VPA2<-function(F.first){
model<-VPA.first(data,M,F.first)
(model$推定漁獲係数[nrow(data),ncol(data)]-model$推定漁獲係数[nrow(data)-1,ncol
(data)])^2
}
#準ニュートン法(quasi-Newton method )による最適化
saiteki<-optim(F.first,VPA2,method="BFGS")
#最適化計算により求められたF.firstを新たに入れ込んで、再計算
VPA.first(data,M,saiteki$par)
}
こんなデータを 使ってみます。これは管理人が適当に作ったものであり、実データではありません。
このデータを使って計算する方法を説明します。Rに詳しい人は以下の説明は読まなくていいです。
このデータがテキストファイル形式でデスクトップ上に保存されていて、data.txt と言う名前だったとします。
Rメニューの ファイル ⇒ ディレクトリの変更 でディレクトリをデスクトップに変更。そしたら
Catch<-read.table("data.txt",header=T)
を実行します。そうしたら、仮想CatchデータがRに読み込まれたことになります。
その後、 このプログラム
をRに張り付けて実行します。（ファイル⇒新しいスクリプト であたらしスクリプトを開いてから張り付けると楽。実行は「F5」キーです。実行は一行ごと
に押さないといけません）そして、
VPA(data=Catch,M=0.2,F.first=0.5)
を実行すると、このよう な結果が出てきます。
VPAにはほかにも色々の計算方法があります。このプログラムはその一例です。
VPAの解説をしてあるPDF資料などはネット上にたくさん落ちています。
このプログラムは、Rらしからぬforループの連続使用により作られています。もっと美しいプログラムならよかったのですが……
