VPAプログラム

★新サイト完成しました!
3秒後に自動的に移動します

変わらない方は こちらからどうぞ
http://logics-of-blue.com/vpa%EF%BC%88%E5%8C%97%E5%A4%A7%E4%BB%95%E6%A7%98%EF%BC%89/

資源量推定手法の王道、VPAのプログラムです。

最終齢はプラスグループで漁業が定常な場合を想定しました。プラスグループの計算においては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ループの連続使用により作られています。もっと美しいプログラムならよかったのですが……


ページトップへ
Copyright (C) 2011- 海と魚と統計解析 All Rights Reserved.