★新サイト完成しました!
３秒後に自動的に移動します
変わらない方は こちらからどうぞ
http://logics-of-blue.com/%E3%83%81%E3%83%A5%E3%83%BC%E3%83%8B%E3%83%B3%E3%82%B0vpa/
資源量推定手法の王道、ＶＰＡのプログラムです。
今回は資源解析手法教科書のpp112－115の チューニングVPAの計算方法をそのまま用いて計算プログラムを作りました。ただし、p115の手法②「6歳のFは（11）式」と書かれていましたが、本 プログラムでは（12）式を使用しました。
VPAの結果を資源量指数でチューニングすること により、解析結果の信頼性を上げようというものですが、資源量指数の正確性に大きく左右されるという欠点もあるようです。細かい解説は教科書に譲ります。
一応教科書p114の表2の数値例を用いて動作確 認を行いましたが、バグ等あればご一報いただけると幸いです。
#===========================
# チューニングVPA プログラム
#===========================
Tuned_VPA<-function(data,index,M,F.first){
data<-data
index<-index
M<-M
F.first<-F.first
#まずは普通に計算
VPA.first<-function(data,M,F.first){
C<-data
#漁獲量
N<-data
#資源尾数
F<-data
#漁獲係数
F[nrow(F)-1,ncol(F)]<-F.first#最高齢一歩手前のFを適当に定める&
nbsp;
M<-M
#自然死亡係数
# 最高齢一歩手前のFから、シングルコホート解析する
N[nrow(N)-1,ncol(N)]<-C[nrow(C)-1,ncol(C)]*( exp(M/2)
)/(1-exp(-1*F[nrow(F)-1,ncol(F)]))
#
↑ まずは最高齢一歩手前魚のNの計算
for (i in
1:(nrow(C)-2)){
#
↓ シングルコホート解析
N[nrow(N)-1-i,ncol(N)-i]<-N[nrow(N)-i,ncol(N)-i+1]*exp(M)+C[nrow(C)-1-i,ncol(C)-i]*exp(M/2)
}
for (i in
1:(nrow(C)-2)){
# ↓ Fの計算
F[nrow(F)-1-i,ncol(F)-i]<-
(-1)*log(1-(
C[nrow(F)-1-i,ncol(F)-i]*exp(M/2) )/N[nrow(F)-1-i,ncol(F)-i])
}
# 次は最高齢Fを近似させてから、横にスライドさせていく
F[nrow(F),ncol(F)]<-F[nrow(F)-1,ncol(F)]
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の近似
}
# 次は近年のFを近似させていく
for (j in (2:(nrow(C)-1))){
F[nrow(F)-j,ncol(F)]<-(F[nrow(F)-j,ncol(F)-1]+F[nrow(F)-j,ncol(F)-2]+F[nrow(F)-j,ncol(F)-3])*
F[nrow(F)-1,ncol(F)]/(F[nrow(F)-1,ncol(F)-1]+F[nrow(F)-1,ncol(F)-2]+F
[nrow(F)-1,ncol(F)-3])
# ↑ 最新年Fを、最新年最高齢一歩手前のFから近似させる
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){
Nsum<-apply(VPA.first(data,M,F.first)$推定資源尾数[,-1],2,sum)
q<-sum(index*Nsum)/sum(Nsum^2)
sum((index-q*Nsum)^2)
}
#最適化
saiteki<-optimise(VPA2,interval=c(0,2))
#最適化計算により求められたF.firstを新たに入れ込んで、再計算
VPA.first(data,M,saiteki$minimum)
}
こ んなデータを 使ってみます。これは管理人が適当に作ったものであり、実データではありません。
また、資源量指数としてはこれを使用しました。かなり雑に作ったサンプルデータです。
このデータを使って計算する方法を説明します。Rに詳しい人は以下の説明は読まなくていいです。
このデータがテキストファイル形式でデスクトップ上に保存されていて、data.txt やindex.txtと言う名前だったとします。
Rメニューの ファイル ⇒ ディレクトリの変更 でディレクトリをデスクトップに変更。そしたら
Catch<-read.table("data.txt",header=T)
Index<-read.table("index.txt",header=T)
を実行します。そうしたら、仮想CatchデータがRに読み込まれたことになります。
その後、 このプログラム
をRに張り付けて実行します。（ファイル⇒新しいスクリプト であたらしスクリプトを開いてから張り付けると楽。実行は「F5」キーです。実行は一行ごと
に押さないといけません）そして、
Tuned_VPA(data=Catch,index=Index,M=0.4,F.first=0.1)
を実行すると、このよう な結果が出てきます。
チューニングVPAにはほかにも色々の計算方法があります。このプログラムはその一例です。
いままでのVPAとは恐ろしく異なる結果となっ
てしまいましたが、資源量指数の作り方がおかしかったからだと思います。 F.firstの値がどんどん0へと近づいて行ってしまったので。実際に残差を
見ると、やはりF.firstが小さい方が残差が小さくなっていました。そのうちシミュレーションをして確かめられれば、と思います。
