knitr::opts_chunk$set( collapse = TRUE, comment = "#>" )
library(frasyr) # 古い関数を呼び出して上書き source("https://raw.githubusercontent.com/ichimomo/future-rvpa/master/future1.11.r") ## 年齢別体重など、将来で仮定する生物パラメータ ## を平均する期間が共通の場合、あらかじめbyearに入れておく ## 同様に、将来で仮定するRPSをとる期間もrps.yearとしておく byear <- 1998:2000 rps.year <- 1991:2000 ## 管理基準値の計算 data(vout3) rout <- ref.F(vout3, waa.year=byear,maa.year=byear, M.year=byear,rps.year=rps.year, max.age=Inf,pSPR=c(20,25,30,40), Fspr.init=0.2) rout$summary # 結果の概観
## 将来予測(vout1); 例えば親魚資源を2007年に400までに回復させる fout1 <- future.vpa(vout3,currentF=NULL, multi=1, nyear=15,start.year=2001,N=1000,ABC.year=2002, waa.year=byear,maa.year=byear,M.year=byear, is.plot=FALSE, rec.new=NULL, recfunc=RPS.simple.rec, rec.arg=list(rps.year=rps.year, upper.ssb=Inf,bias.correction=TRUE,rpsmean=FALSE, upper.recruit=Inf), # ↓回復シナリオ(Frec)の場合、Frecにリストを与えて、回復シナリオを設定する # stochasitc=FALSEで、決定論的な親魚資源量がBlimitと一致するようにする # stochasitc=TRUEでは、確率論的に将来予測をしたとき、50%の確率で # 親魚資源がBlimitを上回るようにする Frec=list(stochastic=FALSE,future.year=2007, Blimit=500000,seed=1,method="optim",scenario="blimit")) fout1$ABC[1] #1000回シミュレーションで1001個の結果が返される。1個目の要素は、決定論的 #将来予測の結果。2-1001個目の要素が、確率論的な将来予測の結果 hist(fout1$ABC[-1]) # 加入の不確実性を考慮した場合のABCの分布 fout1$multi # 最終的に使われたFcurrentに対する乗数 ## 将来予測(vout3);Fmedで漁獲 par(mfrow=c(2,2)) fout2 <- future.vpa(vout3, # ↓NULLの場合、vout1$Fc.at.aが用いられる currentF=NULL, # ↓管理がスタートする年からcurrentFに乗じられる係数 # ここでは,Fmedを使っている multi=rout$summary$Fmed[3], is.plot=FALSE, nyear=15,start.year=2001,N=1000, # ↓ABCを計算する年 ABC.year=2002, waa.year=byear,maa.year=byear,M.year=byear, rec.new=NULL, # ↓将来の加入関数とその引数 recfunc=RPS.simple.rec, rec.arg=list(rps.year=rps.year, upper.ssb=Inf,bias.correction=TRUE,rpsmean=FALSE, upper.recruit=Inf)) fout2$ABC[1] hist(fout2$ABC[-1]) ## ABC計算 SSBcur.sim <- rev(colSums(vout3$ssb))[1] ABC.sim <- getABC(res.vpa=vout3, # vpaの結果 res.ref=rout, # ref.Fの結果 res.future=fout2, # future.vpaの結果 target.year=2005, # 確率を計算する基準の年 is.plot=FALSE, Blim=200,N=1000,SSBcur=SSBcur.sim, # ↓ABCの基礎となる管理基準値(names(rout2$summary)から選ぶ) ref.case=c("Fcurrent","Fcurrent","FpSPR.20.SPR","FpSPR.30.SPR"), multi=c(1,fout2$multi,1,1)) # 上の管理基準値に対する乗数 ABC.sim$ABC # 要約表 ## 結果をcsvファイル、グラフとして出力する out.vpa(res=vout3, # VPAの結果 rres=rout, # 管理基準値の計算結果 fres=fout1, # 将来予測結果 ABC=ABC.sim) # ABC計算結果 # デフォルトではvpa.pdfとvpa.csvというファイルが作成され、 # そこに推定された資源量などの情報が載っています # res, rres, fres, ABCがすべて揃っていなくてもOK # ない場合はNULLを入れる # 例:VPAの結果だけ出力する場合 out.vpa(res=vout3)
## bootstrap結果を使って将来予測 fssb <- fbiom <- ABC <- NULL tmp <- fout2$input tmp$N <- 20 # stochastic runは20回づつ繰り返す for(i in 1:length(boot.sim1)){ tmp$res0 <- boot.sim1[[i]] tmp$res0$input <- vout3$input # ↓ブートストラップで推定に失敗している場合、エラーが出て止まるのでtryで囲む ftmp <- try(do.call(future.vpa,tmp)) if(class(ftmp)!="try-error"){ fssb <- cbind(fssb,ftmp$vssb[,-1]) fbiom <- cbind(fbiom,ftmp$vbiom[,-1]) ABC <- cbind(ABC,ftmp$ABC) } } par(mfrow=c(2,1)) boxplot(t(fssb),ylim=c(0,500)) boxplot(t(fbiom))
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.