★新サイト完成しました!
３秒後に自動的に移動します
変わらない方は こちらからどうぞ
http://logics-of-blue.com/vpa%E8%B3%87%E6%BA%90%E8%A7%A3%E6%9E%90%E6%89%8B%E6%B3%95%E6%95%99%E7%A7%91%E6%9B%B8/
資源量推定手法の王道、ＶＰＡのプログラムです。
前々から載せていたものは、VPAとしては間違っ てはいないのですが、資源解析手法教科書に載っていたものとはややことなる計算式を使用していましたので、ここでは、資源解析手法教科書に載っている計算 式をそのまま用いたプログラムを載せました。
具体的にいうと、資源解析手法教科書
pp110‐112の4節「実用的なVPA」を参考としました。
一応、教科書の問題3の数値例（p107）と表2（p114）の数値を用いて動作確認を行いましたが、バグ等ございましたらご一報いただけると幸いで
す。
#===========================
# VPA プログラム
#===========================
VPA.2<-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]*(
exp(M/2) )/(1-exp(-1*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]<-
(-1)*log(1-( C[nrow(F)-i,ncol(F)-i-j]*exp(M/2)
)/N[nrow(F)-i,ncol(F)-i-j])
}
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)]*(
exp(M/2) )/(1-exp(-1*F[nrow(F)-j,ncol(F)]))
# ↑ 近似した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]<-(-1)*log(1-(
C[nrow(F)-i-j,ncol(F)-i]*exp(M/2) )/N[nrow(F)-i-j,ncol(F)-i])
}
}
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.2(data=Catch,M=0.2,F.first=0.5)
を実行すると、このよう な結果が出てきます。
余談ですが、Mの値を0.8とかとんでもない値を突っ込んでみると、資源量が負の値になってしまいました。エラーとしてそのような値を 表示させないようにしたかったのですが。
こ れは、最適化計算の結果Fが負の値になってしまい、その結果1-exp(-F)の値がマイナスになって、最高齢魚の資源尾数が負になることが原因だと思い ます。 サンプルデータが悪かったようで、教科書のデータであれば、Mの値を多少変えても負にはなりませんでした（M=２とかひどい（ありえない）値を入 れれば負になった）。
