#---------------- 管理基準値計算のための関数 ------------------------
# ref.F
ref.F <- function(
res, # VPAの結果のオブジェクト
sel=NULL, # 仮定する選択率.NULLの場合,res$Fc.at.ageが使われる
waa=NULL, # 仮定する生物パラメータ.直接の値を入れるか,年を指定するやり方のどちらでも動く。直接指定するほうが優先。
maa=NULL,
M=NULL,
waa.catch=NULL,
M.year=NULL,
waa.year=NULL, # 年を指定して生物パラメータを仮定する場合.年の範囲の平均値が用いられる.NULLの場合,VPA最終年の値が使われる
maa.year=NULL,
rps.year = NULL, # Fmedの計算に使うRPSの年の範囲.NULLの場合,全範囲が用いられる
max.age = Inf, # 加入年齢を0歳としたときに、SPR計算で考慮される最大の年齢(年齢の数ではないことに注意)。加入年齢が1歳以上のときは、SPR計算で考慮したい年齢-加入年齢を入力する、またはmin.ageの引数に加入年齢を設定する。
min.age = 0, # 加入年齢が0歳でないときに指定できる
d = 0.001,
Fspr.init = 0.5, # F%SPRの初期値
Fmax.init = 1.5, # Fmaxの初期値
F0.1.init = 0.7, # F0.1の初期値
pSPR = seq(10,90,by=10), # F%SPRを計算するときの%SPR
iterlim=1000,
plot=TRUE,
Pope=FALSE, # 2014.7.4追加
F.range = seq(from=0,to=2,length=101) # YPR, SPR曲線を書くときのFの範囲
){
argname <- ls()
arglist <- lapply(argname,function(x) eval(parse(text=x)))
names(arglist) <- argname
naa <- res$naa
ssb <- res$ssb
ny <- ncol(naa)
years <- dimnames(naa)[[2]]
ages <- dimnames(naa)[[1]]
if(is.null(sel)){
Fc.at.age <- res$Fc.at.age
sel <- Fc.at.age/max(Fc.at.age,na.rm=TRUE)
}
else{
Fc.at.age <- sel
}
sel <- sel/max(sel,na.rm=T)
na <- sum(!is.na(sel))
if(is.null(waa.year)) waa.year <- rev(years)[1]
if(is.null(maa.year)) maa.year <- rev(years)[1]
if(is.null(M.year)) M.year <- rev(years)[1]
if(is.null(rps.year)) rps.year <- as.numeric(colnames(res$naa))
if(is.null(waa)) waa <- apply(as.matrix(as.data.frame(res$input$dat$waa)[as.character(waa.year)]),1,mean)
if(is.null(M)) M <- apply(as.matrix(as.data.frame(res$input$dat$M)[as.character(M.year)]),1,mean)
if(is.null(maa)) maa <- apply(as.matrix(as.data.frame(res$input$dat$maa)[as.character(maa.year)]),1,mean)
if(is.null(waa.catch)){
if(is.null(res$input$dat$waa.catch)){
waa.catch <- waa
}
else{
waa.catch <- apply(as.matrix(as.data.frame(res$input$dat$waa.catch)[as.character(waa.year)]),1,mean)
}
}
ssb.coef <- ifelse(is.null(res$ssb.coef),0,res$ssb.coef)
min.age <- min(as.numeric(rownames(res$naa)))
if(min.age==0) slide.tmp <- TRUE else slide.tmp <- -1:-min.age
rps.data <- data.frame(year=as.numeric(names(colSums(ssb,na.rm=T))),
ssb=as.numeric(colSums(ssb,na.rm=T)),
recruit=as.numeric(c(naa[1,slide.tmp],rep(NA,min.age))))
if (sum(is.na(rps.data$year))>0) rps.data <- rps.data[-which(is.na(rps.data$year)),]
rps.data$rps <- rps <- rps.data$recruit/rps.data$ssb
# rps <- as.numeric(naa[1,]/colSums(ssb,na.rm=TRUE))
# if (is.null(rps.year)) rps.year <- years
tmp <- rps.data$year %in% rps.year
rps.q <- quantile(rps[tmp], na.rm=TRUE, probs=c(0.1,0.5,0.9))
rps.q <- c(rps.q,mean(rps[tmp], na.rm=TRUE))
names(rps.q)[4] <- "mean"
spr.q <- 1/rps.q
# browser()
# F.spr
spr.f.est <- function(log.p, out=FALSE, sub="med", spr0=NULL){
Fr <- exp(log.p)
tmp <- calc.rel.abund(sel,Fr,na,M,waa,waa.catch,maa,min.age=min.age,max.age=max.age,Pope=Pope,ssb.coef=ssb.coef)
rel.abund <- tmp$rel.abund
spr <- sum(tmp$spr)
if (isTRUE(out)) obj <- spr
else{
if(sub=="mean") obj <- (spr-spr.q[4])^2
if(sub=="low") obj <- (spr-spr.q[3])^2
if(sub=="med") obj <- (spr-spr.q[2])^2
if(sub=="high") obj <- (spr-spr.q[1])^2
if(is.numeric(sub)) obj <- (spr/spr0-sub/100)^2
}
return(obj)
}
spr0 <- spr.f.est(-Inf, out=TRUE)
Fmed.res <- nlm(spr.f.est, Fspr.init, out=FALSE, sub="med", iterlim = iterlim)
Fmean.res <- nlm(spr.f.est, Fspr.init, out=FALSE, sub="mean", iterlim = iterlim)
# browser()
Flow.res <- nlm(spr.f.est, Fspr.init, out=FALSE, sub="low", iterlim = iterlim)
Fhigh.res <- nlm(spr.f.est, Fspr.init, out=FALSE, sub="high", iterlim = iterlim)
Fmean <- exp(Fmean.res$estimate)
Fmed <- exp(Fmed.res$estimate)
Flow <- exp(Flow.res$estimate)
Fhigh <- exp(Fhigh.res$estimate)
if (!is.null(pSPR)){
FpSPR <- NULL
for (i in pSPR){
FpSPR.res <- nlm(spr.f.est, Fspr.init, out=FALSE, sub=i, spr0=spr0, iterlim=iterlim)
# print(FpSPR.res)
# cat("i", FpSPR.res$code," ")
FpSPR <- c(FpSPR, exp(FpSPR.res$estimate))
}
names(FpSPR) <- paste(pSPR,"%SPR",sep="")
}
# Fmax
ypr.f.est <- function(log.p, out=FALSE){
Fr <- exp(log.p)
tmp <- calc.rel.abund(sel,Fr,na,M,waa,waa.catch,maa,max.age=max.age,Pope=Pope,ssb.coef=ssb.coef)
rel.abund <- tmp$rel.abund
ypr <- sum(tmp$ypr)
if (isTRUE(out)) obj <- ypr else obj <- -ypr
return(obj)
}
Fmax.res <- nlm(ypr.f.est, log(Fmax.init), out=FALSE)
Fmax <- exp(Fmax.res$estimate)
# F0.1
Fp <- function(log.p, out=FALSE){
Fr <- exp(log.p)
tmp <- calc.rel.abund(sel,Fr,na,M,waa,waa.catch,maa,max.age=max.age,Pope=Pope,ssb.coef=ssb.coef)
rel.abund <- tmp$rel.abund
ypr <- sum(tmp$ypr)
if (isTRUE(out)) obj <- ypr else obj <- -ypr
return(obj)
}
F0.1.est <- function(log.p){
p <- exp(log.p)
ref.trend <- (Fp(log(d))-Fp(log(0)))/d
trend <- (Fp(log(p+d)) - Fp(log(p)))/d
obj <- (ref.trend/10 - trend)^2
obj
}
F0.1.res <- nlm(F0.1.est,log(F0.1.init))
F0.1 <- exp(F0.1.res$estimate)
# Fcurrent
Fcurrent <- c(max(Fc.at.age,na.rm=T), mean(Fc.at.age,na.rm=T))
# output
f.mean <- function(x) mean(x*sel, na.rm=T)
Fmean <- c(Fmean, f.mean(Fmean))
Fmed <- c(Fmed, f.mean(Fmed))
Flow <- c(Flow, f.mean(Flow))
Fhigh <- c(Fhigh, f.mean(Fhigh))
Fmax <- c(Fmax, f.mean(Fmax))
F0.1 <- c(F0.1, f.mean(F0.1))
names(Fcurrent) <- names(Fmed) <- names(Fmean) <- names(Flow) <- names(Fhigh) <- names(Fmax) <- names(F0.1) <- c("max","mean")
Res <- list(sel=sel, min.age=min.age, max.age=max.age, rps.q=rps.q, spr.q=spr.q, Fcurrent=Fcurrent, Fmed=Fmed, Flow=Flow, Fhigh=Fhigh, Fmax=Fmax, F0.1=F0.1, Fmean=Fmean,rps.data=rps.data)
if (!is.null(pSPR)){
FpSPR <- rbind(FpSPR, sapply(FpSPR, f.mean))
rownames(FpSPR) <- c("max","mean")
Res$FpSPR <- FpSPR
}
#---- make summary
Res$summary <- as.data.frame(Res[substr(names(Res),1,1)=="F"])
Res$summary <- rbind(Res$summary,Res$summary[1,]/Res$summary[1,1])
dimnames(Res$summary)[[1]][3] <- "Fref/Fcur"
#----- YPR & SPR figure -----
F_current <- Res$summary$Fcurrent[1]*Res$summary$Fcurrent[3]
F.range <- sort(c(F.range, F_current))
spr0 <- sum(calc.rel.abund(sel,0,na,M,waa,waa.catch,maa,min.age=min.age,max.age=max.age,Pope=Pope,ssb.coef=ssb.coef)$spr)
tmp <- lapply(F.range, function(x) calc.rel.abund(sel,x,na,M,waa,waa.catch,maa,min.age=min.age,max.age=max.age,Pope=Pope,ssb.coef=ssb.coef))
ypr <- sapply(tmp,function(x) sum(x$ypr))
spr <- sapply(tmp,function(x) sum(x$spr))/spr0*100
Res$ypr.spr <- data.frame(F.range=F.range,ypr=ypr,spr=spr)
Res$waa <- waa
Res$waa.catch <- waa.catch
Res$maa <- maa
#------------------------------
Res$arglist <- arglist
Res$spr0 <- spr0
Res$ypr.spr$Frange2Fcurrent <- Res$ypr.spr$F.range/F_current
class(Res) <- "ref"
if(isTRUE(plot)){
plot.Fref(Res)
}
return(Res)
}
plot.Fref <- function(rres,xlabel="max", # or, "mean","Fref/Fcur"
vline.text=c("F0.1","Fmax","Fcurrent","Fmed") # and "FpSPR.20.SPR" etc..
){
old.par <- par()
par(mar=c(4,4,1,4))
F.range <- rres$ypr.spr$F.range
if(xlabel=="Fref/Fcur") F.range <- F.range/rres.pma$summary$Fcurrent[1]*rres.pma$summary$Fcurrent[3]
if(xlabel=="mean") F.range <- F.range/rres.pma$summary$Fcurrent[1]*rres.pma$summary$Fcurrent[2]
spr <- rres$ypr.spr$spr
ypr <- rres$ypr.spr$ypr
plot(F.range,spr,xlab=xlabel,ylab="%SPR",type="l",ylim=c(0,max(spr)))
par(new=T)
plot(F.range,ypr,axes=F,xlab="",ylab="",lty=2,type="l",ylim=c(0,max(ypr)))
axis(side=4)
mtext("YPR",side=4,line=2)
n.line <- which(rownames(rres$summary) %in% xlabel)
abline(v=xx <- c(rres$summary[vline.text][n.line,]))
text(xx,max(ypr)*seq(from=0.5,to=0.3,length=length(vline.text)),vline.text)
legend("topright",lty=1:2,legend=c("SPR","YPR"))
options(warn=-1); par(old.par); options(warn=0)
}
calc.rel.abund <- function(sel,Fr,na,M,waa,waa.catch=NULL,maa,min.age=0,max.age=Inf,Pope=TRUE,ssb.coef=0){
if(is.null(waa.catch)) waa.catch <- waa
rel.abund <- rep(NA, na)
rel.abund[1] <- 1
for (i in 2:(na-1)) {
rel.abund[i] <- rel.abund[i-1]*exp(-M[i-1]-sel[i-1]*Fr)
}
rel.abund[na] <- rel.abund[na-1]*exp(-M[na-1]-sel[na-1]*Fr)*(1-exp(-((max.age-min.age)-(na-2))*(M[na]+sel[na]*Fr)))/(1-exp(-M[na]-sel[na]*Fr))
if(isTRUE(Pope)){
ypr1 <- rel.abund*waa.catch[1:na]*(1-exp(-sel[1:na]*Fr))*exp(-M[1:na]/2)
}
else{
# use Baranov catch equation
ypr1 <- rel.abund*(1-exp(-sel[1:na]*Fr-M[1:na]))*sel[1:na]*Fr/
(sel[1:na]*Fr+M[1:na])*waa.catch[1:na]
}
spr <- rel.abund*waa[1:na]*maa[1:na]*exp(-ssb.coef*(sel[1:na]*Fr+M[1:na]))
return(list(rel.abund=rel.abund,ypr=ypr1,spr=spr))
}
##----------------------- 将来予測関数 ----------------------------
## multiのオプションは管理後のFのmultiplier(管理前後でselectivityが同じ)
future.vpa <-
function(res0,
currentF=NULL, # 管理前のF
multi=1, # 管理後(ABC.yearから)のF (current F x multi)
nyear=10,Pope=res0$input$Pope,
outtype="FULL",
multi.year=1,#ある特定の年だけFを変えたい場合。デフォルトは1。変える場合は、指定した年またはタイムステップの要素数のベクトルで指定。
# 年数の指定
start.year=NULL, # 将来予測の開始年,NULLの場合はVPA計算の最終年の次の年
ABC.year=NULL, # ABC yearを計算する年。NULLの場合はVPA計算の最終年の次の次の年
waa.year=NULL, # VPA結果から生物パラメータをもってきて平均する期間
# NULLの場合,VPAの最終年のパラメータを持ってくる
maa.year=NULL, # VPA結果から生物パラメータをもってきて平均する期間
M.year=NULL, # VPA結果から生物パラメータをもってきて平均する期間
seed=NULL,
strategy="F", # F: 漁獲係数一定, E: 漁獲割合一定、C: 漁獲量一定(pre.catchで漁獲量を指定)
HCR=NULL,# HCRを使う場合、list(Blim=154500, Bban=49400,beta=1,year.lag=0)のように指定するか、以下の引数をセットする,year.lag=0で将来予測年の予測SSBを使う。-2の場合は2年遅れのSSBを使う
use.MSE=FALSE,
beta=NULL,delta=NULL,Blim=0,Bban=0,
plus.group=res0$input$plus.group,
N=1000,# 確率的なシミュレーションをする場合の繰り返し回数。
# N+1の結果が返され、1列目に決定論的な結果が
# 0を与えると決定論的な結果のみを出力
silent=FALSE, is.plot=TRUE, # 計算条件を出力、プロットするか
random.select=NULL, # 選択率をランダムリサンプリングする場合、ランダムリサンプリングする年を入れる
# strategy="C"または"E"のときのみ有効
pre.catch=NULL, # list(year=2012,wcatch=13000), 漁獲重量をgivenで与える場合
# list(year=2012:2017,E=rep(0.5,6)), 漁獲割合をgivenで与える場合
##-------- 加入に関する設定 -----------------
rec.new=NULL, # 指定した年の加入量
# 年を指定しないで与える場合は、自動的にスタート年の加入になる。
# list(year=, rec=)で与える場合は、対応する年の加入を置き換える。
##--- 加入関数
recfunc=HS.recAR, # 再生産関係の関数
rec.arg=list(a=1,b=1,rho=0,sd=0,c=1,bias.correction=TRUE,
resample=FALSE,resid=0,resid.year=NULL), # 加入の各種設定
##--- Frecオプション;Frec計算のための設定リストを与えると、指定された設定でのFrecに対応するFで将来予測を行う
Frec=NULL,
# list(stochastic=TRUE, # TRUEの場合、stochastic simulationで50%の確率でBlimitを越す(PMS, TMI)
# FALSEの場合、RPS固定のprojectionがBilmitと一致する(NSK)
# future.year=2018, # 何年の資源量を見るか?
# Blimit=450*1000, # Blimit (xトン)
# scenario="catch.mean" or "blimit" (デフォルトはblimit; "catch.mean"とするとstochastic simulationにおける平均漁獲量がBlimitで指定した値と一致するようになる)
# Frange=c(0.01,2*mult)) # Fの探索範囲
waa=NULL,waa.catch=NULL,maa=NULL,M=NULL, # 季節毎の生物パラメータ、または、生物パラメータを外から与える場合
replace.rec.year=2012, # 加入量を暦年の将来予測での加入量に置き換えるか?
F.sigma=0,
waa.fun=FALSE, #waaをnaaのfunctionとするか
naa0=NULL,eaa0=NULL,ssb0=NULL,faa0=NULL,
add.year=0, # 岡村オプションに対応。=1で1年分余計に計算する
det.run=TRUE # 1回めのランは決定論的将来予測をする(完璧には対応していない)
){
argname <- ls()
arglist <- lapply(argname,function(x) eval(parse(text=x)))
names(arglist) <- argname
if(is.null(res0$input$unit.waa)) res0$input$unit.waa <- 1
if(is.null(res0$input$unit.caa)) res0$input$unit.caa <- 1
if(is.null(res0$input$unit.biom)) res0$input$unit.biom <- 1
if(is.null(plus.group)) plus.group <- TRUE
if(is.null(Pope)) Pope <- FALSE
##--------------------------------------------------
if(isTRUE(det.run)) N <- N + 1
years <- as.numeric(dimnames(res0$naa)[[2]])
##------------- set default options
if(is.null(currentF)) currentF <- res0$Fc.at.age
if(is.null(waa.year)) waa.year <- rev(years)[1]
if(is.null(maa.year)) maa.year <- rev(years)[1]
if(is.null(M.year)) M.year <- rev(years)[1]
if(is.null(start.year)) start.year <- rev(years)[1]+1
if(is.null(ABC.year)) ABC.year <- rev(years)[1]+1
## if(!is.null(Bban)) Bban$is.Bban <- rep(FALSE,N)
arglist$ABC.year <- ABC.year
##-------------
##---- set S-R functin option -----
## 使う関数によっては必要ないオプションもあるが、使わないオプションを入れてもエラーは出ないので、
# rec.arg$resampleがNULLかどうかで、パラメトリックな誤差分布かそうでないか(残差リサンプリング)を判別する
if(is.null(rec.arg$rho)){
rec.arg$rho <- 0
if(!silent) cat("rec.arg$rho is assumed to be 0...\n")
}
if(is.null(rec.arg$sd2)) rec.arg$sd2 <- sqrt(rec.arg$sd^2/(1-rec.arg$rho^2)) #rho込み平均補正用SD # HS.recAR
## resampling optionを使わない場合
if(is.null(rec.arg$resample)|!isTRUE(rec.arg$resample)){
if(is.null(rec.arg$bias.correction)) rec.arg$bias.correction <- TRUE # HS.recAR, HS.rec0
if(is.null(rec.arg$rho)){
rec.arg$rho <- 0 # HS.recAR, HS.rec0
rec.arg$resid <- 0
}
if(!is.null(rec.arg$rho)){
if(rec.arg$rho>0){
if(is.null(eaa0)){
if(is.null(rec.arg$resid.year)) rec.arg$resid <- rep(rev(rec.arg$resid)[1],N)
else rec.arg$resid <- rep(mean(rev(rec.arg$resid)[1:rec.arg$resid.year]),N)
}
else{
rec.arg$resid <- eaa0
}
}
else{
rec.arg$resid <- rep(0,N)
}
}
}
else{
if(rec.arg$rho>0) stop("You set rho is >0. You cannot use resample=TRUE option when rho>0") # resamplingの場合に自己相関は考慮できないのでrhoは強制的にゼロ
}
if(!is.null(rec.arg$sd)) rec.arg$sd <- c(0,rep(rec.arg$sd,N-1))
if(!is.null(rec.arg$sd2)) rec.arg$sd2 <- c(0,rep(rec.arg$sd2,N-1))
if(!is.null(HCR) && is.null(HCR$year.lag)) HCR$year.lag <- 0
##---------------------------------
if(!is.null(beta)){
HCR$beta <- beta
HCR$Blim <- Blim
HCR$Bban <- Bban
}
# fyears <- seq(from=start.year,to=start.year+nyear-1,by=1/ts)
fyears <- seq(from=start.year,to=start.year+nyear+add.year,by=1)
fyear.year <- floor(fyears)
ntime <- length(fyears)
ages <- as.numeric(dimnames(res0$naa)[[1]]) # ages:VPAで考慮される最大年齢数
min.age <- min(as.numeric(ages))
year.overlap <- years %in% start.year
{if(sum(year.overlap)==0){
nage <- sum(!is.na(res0$naa[,ncol(res0$naa)])) # nage:将来予測で考慮すべき年の数
}
else{
nage <- sum(!is.na(res0$naa[,year.overlap]))
}}
if(!silent){
arglist.tmp <- arglist
arglist.tmp$res0 <- NULL
arglist.tmp$Bban <- arglist.tmp$Bblim <- arglist.tmp$beta <- arglist.tmp$ssb0 <- arglist.tmp$strategy <- NULL
print(arglist.tmp)
}
# シードの設定
if(is.null(seed)) arglist$seed <- as.numeric(Sys.time())
#------------Frecオプションの場合 -------------
if(!is.null(Frec)){
multi.org <- multi
if(is.null(Frec$stochastic)) Frec$stochastice <- TRUE
if(is.null(Frec$target.probs)) Frec$target.probs <- 50
if(is.null(Frec$scenario)) Frec$scenario <- "blimit" # 2017/12/25追記
if(is.null(Frec$Frange)) Frec$Frange <- c(0.01,multi.org*2) # 2017/12/25追記(探索するFの範囲の指定)
if(is.null(Frec$future.year)) Frec$future.year <- fyears[length(fyears)]-1
# arglist$Frec <- Frec
getFrec <- function(x,arglist){
set.seed(arglist$seed)
arglist.tmp <- arglist
arglist.tmp$multi <- x
arglist.tmp$silent <- TRUE
arglist.tmp$Frec <- NULL
arglist.tmp$is.plot <- FALSE
if(Frec$stochastic==FALSE){
arglist.tmp$N <- 0
}
fres.tmp <- do.call(future.vpa,arglist.tmp)
tmp <- rownames(fres.tmp$vssb)==Frec$future.year
if(all(tmp==FALSE)) stop("nyear should be longer than Frec$future.year.")
if(Frec$stochastic==TRUE){
if(Frec$scenario=="blimit"){
is.lower.ssb <- fres.tmp$vssb<Frec$Blimit
probs <- (sum(is.lower.ssb[tmp,-1],na.rm=T)-1)/
(length(is.lower.ssb[tmp,-1])-1)*100
return.obj <- probs-Frec$target.probs
}
# stochastic projectionにおける平均漁獲量を目的の値に一致させる
if(Frec$scenario=="catch.mean"){
return.obj <- (log(Frec$Blimit)-log(mean(fres.tmp$vwcaa[tmp,-1])))^2
}
# stochastic projectionにおける平均親魚資源量を目的の値に一致させる
if(Frec$scenario=="ssb.mean"){
return.obj <- (log(Frec$Blimit)-log(mean(fres.tmp$vssb[tmp,-1])))^2
}
}
else{
return.obj <- Frec$Blimit-fres.tmp$vssb[tmp,1]
}
# return(ifelse(Frec$method=="nibun",return.obj,return.obj^2))
return(return.obj^2)
}
res <- optimize(getFrec,interval=Frec$Frange,arglist=arglist)
multi <- res$minimum
cat("F multiplier=",multi,"\n")
}
#-------------- main function ---------------------
waa.org <- waa
waa.catch.org <- waa.catch
maa.org <- maa
M.org <- M
if(strategy=="C"|strategy=="E") multi.catch <- multi else multi.catch <- 1
faa <- naa <- waa <- waa.catch <- maa <- M <- caa <-
array(NA,dim=c(length(ages),ntime,N),dimnames=list(age=ages,year=fyears,nsim=1:N))
allyears <- sort(unique(c(fyears,years)))
# 全部のデータを記録したフレーム
naa_all <- waa_all <- waa_catch_all <- maa_all <- array(NA,dim=c(length(ages),length(allyears),N),dimnames=list(age=ages,year=allyears,nsim=1:N))
naa_all[,1:length(years),] <- unlist(res0$naa)
waa_all[,1:length(years),] <- unlist(res0$input$dat$waa)
if(is.null(res0$input$dat$waa.catch)) waa_catch_all[,1:length(years),] <- unlist(res0$input$dat$waa) else waa_catch_all[,1:length(years),] <- unlist(res0$input$dat$waa.catch)
maa_all[,1:length(years),] <- unlist(res0$input$dat$maa)
i_all <- which(allyears%in%start.year)
alpha <- thisyear.ssb <- array(1,dim=c(ntime,N),dimnames=list(year=fyears,nsim=1:N))
# future biological patameter
if(!is.null(M.org)) M[] <- M.org else M[] <- apply(as.matrix(res0$input$dat$M[,years %in% M.year]),1,mean)
if(!is.null(waa.org)) waa[] <- waa.org else waa[] <- apply(as.matrix(res0$input$dat$waa[,years %in% waa.year]),1,mean)
if(!is.null(maa.org)) maa[] <- maa.org else maa[] <- apply(as.matrix(res0$input$dat$maa[,years %in% maa.year]),1,mean)
if(!is.null(waa.catch.org)) waa.catch[] <- waa.catch.org
else{
if(!is.null(res0$input$dat$waa.catch)) waa.catch[] <- apply(as.matrix(res0$input$dat$waa.catch[,years %in% waa.year]),1,mean)
else waa.catch <- waa
}
# future F matrix
faa[] <- currentF*multi # *exp(rnorm(length(faa),0,F.sigma))
# ABCyear以前はcurrent Fを使う。
faa[,fyears<min(ABC.year),] <- currentF*exp(rnorm(length(faa[,fyears<min(ABC.year),]),0,F.sigma))
## VPA期間と将来予測期間が被っている場合、VPA期間のFはVPAの結果を使う
overlapped.years <- list(future=which(fyear.year %in% years),vpa=which(years %in% fyear.year))
if(length(overlapped.years$future)>0){
# for(jj in 1:length(vpayears.overlapped)){
for(j in 1:length(overlapped.years$future)){
if(any(res0$faa[,overlapped.years$vpa[j]]>0) && !is.null(res0$input$dat$waa[,overlapped.years$vpa[j]])){ # もしfaaがゼロでないなら(PMIの場合、2012までデータが入っているが、faaはゼロになっているので
faa[,overlapped.years$future[j],] <- res0$faa[,overlapped.years$vpa[j]]
waa[,overlapped.years$future[j],] <- res0$input$dat$waa[,overlapped.years$vpa[j]]
if(!is.null(res0$input$dat$waa.catch)){
waa.catch[,overlapped.years$future[j],] <- res0$input$dat$waa.catch[,overlapped.years$vpa[j]]
}
else{
waa.catch[,overlapped.years$future[j],] <- res0$input$dat$waa[,overlapped.years$vpa[j]]
}
}
}}
#}
tmp <- aperm(faa,c(2,1,3))
tmp <- tmp*multi.year
faa <- aperm(tmp,c(2,1,3))
# vpa.multi <- ifelse(is.null(vpa.mode),1,vpa.mode$multi)
# rps assumption
rps.mat <- array(NA,dim=c(ntime,N),dimnames=list(fyears,1:N))
eaa <- matrix(0,ntime,N)
rec.tmp <- list(rec.resample=NULL,tmparg=NULL)
if (waa.fun){ #年齢別体重の予測関数
WAA <- res0$input$dat$waa
NAA <- res0$naa
# nage <- nrow(WAA)
WAA.res <- lapply(1:nage, function(i) {
log.w <- as.numeric(log(WAA[i,]))
log.n <- as.numeric(log(NAA[i,]))
lm(log.w~log.n)
})
WAA.cv <- sapply(1:nage, function(i) sqrt(mean(WAA.res[[i]]$residuals^2)))
WAA.b0 <- sapply(1:nage, function(i) as.numeric(WAA.res[[i]]$coef[1]))
WAA.b1 <- sapply(1:nage, function(i) as.numeric(WAA.res[[i]]$coef[2]))
## waa.rand <- array(0,dim=c(al,nyear+1-min.age,N))
set.seed(0)
cv.vec <- rep(WAA.cv,N*ntime)
waa.rand <- array(rnorm(length(cv.vec),-0.5*cv.vec^2,cv.vec),dim=c(nage,ntime,N))
waa.rand[,,1] <- 0
}
set.seed(arglist$seed)
# 将来予測の最初の年の設定;バリエーションがありややこしいのでここで設定される
if(!start.year%in%years){
# VPA結果が2011年までで、将来予測の開始年が2012年の場合
if(start.year==(max(years)+1)){
{if(is.null(res0$input$dat$M)){
M.lastyear <- M.org
}
else{
M.lastyear <- res0$input$dat$M[,length(years)]
}}
# 1年分forwardさせた年齢構成を初期値とする
tmp <- forward.calc.simple(res0$faa[1:nage,length(years)],
res0$naa[1:nage,length(years)],
M.lastyear[1:nage],
plus.group=plus.group)
naa[1:nage,1,] <- naa_all[1:nage,i_all,] <- tmp
if(fyears[1]-min.age < start.year){
thisyear.ssb[1,] <- sum(res0$ssb[,as.character(fyears[1]-min.age)],na.rm=T)
# thisyear.ssb <- rep(thisyear.ssb,N)
}
else{
if(waa.fun){
waa[2:nage,1,] <- t(sapply(2:nage, function(ii) as.numeric(exp(WAA.b0[ii]+WAA.b1[ii]*log(naa[ii,1,])+waa.rand[ii,1,]))))
}
thisyear.ssb[1,] <- colSums(naa[,1,]*waa[,1,]*maa[,1,],na.rm=T)*res0$input$unit.waa/res0$input$unit.biom }
thisyear.ssb[1,] <- thisyear.ssb[1,]+(1e-10)
if(!is.null(ssb0)) thisyear.ssb[1,] <- colSums(ssb0)
rec.tmp <- recfunc(thisyear.ssb[1,],res0,
rec.resample=rec.tmp$rec.resample,
rec.arg=rec.arg)
eaa[1,] <- rec.tmp$rec.resample[1:N]
rec.arg$resid <- rec.tmp$rec.resample # ARオプションに対応
if(!is.null(rec.tmp$rec.arg)) rec.arg <- rec.tmp$rec.arg
naa[1,1,] <- naa_all[1,i_all,] <- rec.tmp$rec
if (waa.fun) {
waa[1,1,] <- as.numeric(exp(WAA.b0[1]+WAA.b1[1]*log(naa[1,1,])+waa.rand[1,1,]))
}
rps.mat[1,] <- naa[1,1,]/thisyear.ssb[1,]
}
else{
stop("ERROR Set appropriate year to start projection\n")
}
}
else{
# VPA期間と将来予測期間が被っている場合にはVPAの結果を初期値として入れる
naa[,1,] <- naa_all[,i_all,] <- res0$naa[,start.year==years]
}
# もし引数naa0が与えられている場合にはそれを用いる
if(!is.null(naa0)){
naa[,1,] <- naa_all[,i_all,] <- naa0
if(is.null(faa0)) faa0 <- res0$Fc.at.age
faa[] <- faa0*multi
}
if(!is.null(rec.new)){
if(!is.list(rec.new)){
naa[1,1,] <- naa_all[1,i_all,] <- rec.new
}
else{ # rec.newがlistの場合
naa[1,fyears%in%rec.new$year,] <- naa_all[,allyears%in%rec.new$year,] <- rec.new$rec
}}
# 2年目以降の将来予測
for(i in 1:(ntime-1)){
#漁獲量がgivenの場合
if(!is.null(pre.catch) && fyears[i]%in%pre.catch$year){
if(!is.null(pre.catch$wcatch)){
if(fyears[i]<ABC.year){
tmpcatch <- as.numeric(pre.catch$wcatch[pre.catch$year==fyears[i]])
}
else{
tmpcatch <- as.numeric(pre.catch$wcatch[pre.catch$year==fyears[i]]) * multi.catch
}
}
if(!is.null(pre.catch$E)){
biom <- sum(naa[,i,]*waa[,i,]*res0$input$unit.waa/res0$input$unit.biom)
if(fyears[i]<ABC.year){
tmpcatch <- as.numeric(pre.catch$E[pre.catch$year==fyears[i]]) * biom
}
else{
tmpcatch <- as.numeric(pre.catch$E[pre.catch$year==fyears[i]]) * biom * multi.catch
}
}
# 選択率をランダムサンプリングする場合
# if(!is.null(random.select)) saa.tmp <- as.numeric(res0$saa[,colnames(res0$saa)==sample(random.select,1)])
saa.tmp <- sweep(faa[,i,],2,apply(faa[,i,],2,max),FUN="/")
tmp <- lapply(1:dim(naa)[[3]],
function(x) caa.est.mat(naa[,i,x],saa.tmp[,x],
waa.catch[,i,x],M[,i,x],tmpcatch,Pope=Pope))
faa.new <- sweep(saa.tmp,2,sapply(tmp,function(x) x$x),FUN="*")
caa[,i,] <- sapply(tmp,function(x) x$caa)
faa[,i,] <- faa.new
}
else{
faa.new <- NULL
}
## HCRを使う場合(当年の資源量から当年のFを変更する)
if(!is.null(HCR) && fyears[i]>=ABC.year
&& is.null(faa.new)) # <- pre.catchで漁獲量をセットしていない
{
tmp <- i+HCR$year.lag
if(tmp>0){
ssb.tmp <- colSums(naa[,tmp,]*waa[,tmp,]*maa[,tmp,],na.rm=T)*
res0$input$unit.waa/res0$input$unit.biom
}
else{
vpayear <- fyears[i]+HCR$year.lag
ssb.tmp <- sum(res0$ssb[as.character(vpayear)])
}
alpha[i,] <- ifelse(ssb.tmp<HCR$Blim,HCR$beta*(ssb.tmp-HCR$Bban)/(HCR$Blim-HCR$Bban),HCR$beta)
alpha[i,] <- ifelse(alpha[i,]<0,0,alpha[i,])
faa[,i,] <- sweep(faa[,i,],2,alpha[i,],FUN="*")
faa[,i,] <- ifelse(faa[,i,]<0,0,faa[,i,])
}
if(isTRUE(use.MSE)){
aa <- get_ABC_inMSE(naa_all,waa_all,maa_all,faa,M,res0,
i_all-2,2,recfunc,rec.arg,Pope,HCR,plus.group,min.age)
}
## 漁獲して1年分前進(加入はまだいれていない)
tmp <- forward.calc.mat2(faa[,i,],naa[,i,],M[,i,],plus.group=plus.group)
# 既に値が入っているところ(1年目の加入量)は除いて翌年のNAAを入れる
naa.tmp <- naa[,i+1,]
naa.tmp[is.na(naa.tmp)] <- tmp[is.na(naa.tmp)]
naa[,i+1, ] <- naa_all[,i_all+1,] <- naa.tmp
## 当年の加入の計算
if(fyears[i+1]-min.age < start.year){
# 参照する親魚資源量がVPA期間である場合、VPA期間のSSBをとってくる
thisyear.ssb[i+1,] <- sum(res0$ssb[,as.character(fyears[i+1]-min.age)],na.rm=T)*res0$input$unit.waa/res0$input$unit.biom
# thisyear.ssb <- rep(thisyear.ssb,N)
if(!is.null(ssb0)) thisyear.ssb[i+1,] <- colSums(ssb0)
}
else{
# そうでない場合
if(waa.fun){
# 動的なwaaは対応する年のwaaを書き換えた上で使う?
waa[2:nage,i+1-min.age,] <- t(sapply(2:nage, function(ii) as.numeric(exp(WAA.b0[ii]+WAA.b1[ii]*log(naa[ii,i+1-min.age,])+waa.rand[ii,i+1-min.age,]))))
}
thisyear.ssb[i+1,] <- colSums(naa[,i+1-min.age,]*waa[,i+1-min.age,]*maa[,i+1-min.age,],na.rm=T)*res0$input$unit.waa/res0$input$unit.biom
}
thisyear.ssb[i+1,] <- thisyear.ssb[i+1,]+(1e-10)
rec.tmp <- recfunc(thisyear.ssb[i+1,],res0,
rec.resample=rec.tmp$rec.resample,
rec.arg=rec.arg)
if(is.na(naa[1,i+1,1])) naa[1,i+1,] <- naa_all[1,i_all+1,] <- rec.tmp$rec
if (waa.fun) {
waa[1,i+1,] <- as.numeric(exp(WAA.b0[1]+WAA.b1[1]*log(naa[1,i+1,])+waa.rand[1,i+1,]))
}
# if(!is.null(rec.tmp$rec.arg)) rec.arg <- rec.tmp$rec.arg
rps.mat[i+1,] <- naa[1,i+1,]/thisyear.ssb[i+1,]
eaa[i+1,] <- rec.tmp$rec.resample[1:N]
rec.arg$resid <- rec.tmp$rec.resample # ARオプションに対応
i_all <- i_all+1
}
if (!is.null(rec.arg$rho)) rec.tmp$rec.resample <- NULL
if(Pope){
caa[] <- naa*(1-exp(-faa))*exp(-M/2)
}
else{
caa[] <- naa*(1-exp(-faa-M))*faa/(faa+M)
}
caa <- caa[,-ntime,,drop=F]
if(isTRUE(waa.fun)){ ## アドホックな対応! waa.fun=TRUEかつwaa.catchが与えられているとき動かない。また、pre.catchが与えられていてwaa.fun=TRUEの場合も不具合おこる!
waa.catch <- waa[,-ntime,,drop=F]
}
else{
waa.catch <- waa.catch[,-ntime,,drop=F]
}
thisyear.ssb <- thisyear.ssb[-ntime,,drop=F]
waa <- waa[,-ntime,,drop=F]
maa <- maa[,-ntime,,drop=F]
naa <- naa[,-ntime,,drop=F]
faa <- faa[,-ntime,,drop=F]
alpha <- alpha[-ntime,,drop=F]
M <- M[,-ntime,,drop=F]
fyears <- fyears[-ntime]
biom <- naa*waa*res0$input$unit.waa/res0$input$unit.biom
ssb <- naa*waa*maa*res0$input$unit.waa/res0$input$unit.biom
wcaa <- caa*waa.catch*res0$input$unit.waa/res0$input$unit.biom
vwcaa <- apply(wcaa,c(2,3),sum,na.rm=T)
ABC <- apply(as.matrix(vwcaa[fyears%in%ABC.year,,drop=F]),2,sum)
if(!is.null(rec.arg$resample)) if(rec.arg$resample==TRUE) eaa[] <- NA # resamplingする場合にはeaaにはなにも入れない
fres <- list(faa=faa,naa=naa,biom=biom,baa=biom,ssb=ssb,wcaa=wcaa,caa=caa,M=M,rps=rps.mat,
maa=maa,vbiom=apply(biom,c(2,3),sum,na.rm=T),
recruit=naa[1,,],
eaa=eaa,alpha=alpha,thisyear.ssb=thisyear.ssb,
waa=waa,waa.catch=waa.catch,currentF=currentF,
vssb=apply(ssb,c(2,3),sum,na.rm=T),vwcaa=vwcaa,naa_all=naa_all,
years=fyears,fyear.year=fyear.year,ABC=ABC,recfunc=recfunc,rec.arg=rec.arg,
waa.year=waa.year,maa.year=maa.year,multi=multi,multi.year=multi.year,
Frec=Frec,rec.new=rec.new,pre.catch=pre.catch,input=arglist)
if(is.plot){
par(mfrow=c(2,2))
plot.future(fres)
}
if(waa.fun) fres$waa.reg <- WAA.res
if(outtype=="Det"){
fres <- list(faa=faa[,,1],M=M[,,1],recruit=naa[1,,],eaa=eaa,baa=biom,
maa=maa[,,1],vbiom=apply(biom,c(2,3),sum,na.rm=T),
waa=waa[,,1],waa.catch=waa.catch[,,1],currentF=currentF,
vssb=apply(ssb,c(2,3),sum,na.rm=T),vwcaa=vwcaa,alpha=alpha,
years=fyears,fyear.year=fyear.year,ABC=ABC,recfunc=recfunc,
waa.year=waa.year,maa.year=maa.year,multi=multi,multi.year=multi.year,
Frec=Frec,rec.new=rec.new,pre.catch=pre.catch,input=arglist)
}
if(outtype=="short"){
fres <- list(recruit=naa[1,,],eaa=eaa,#baa=biom,
vbiom=apply(biom,c(2,3),sum,na.rm=T),
currentF=currentF,alpha=alpha,
vssb=apply(ssb,c(2,3),sum,na.rm=T),vwcaa=vwcaa,
years=fyears,fyear.year=fyear.year,ABC=ABC,
waa.year=waa.year,maa.year=maa.year,multi=multi,multi.year=multi.year,
Frec=Frec,rec.new=rec.new,pre.catch=pre.catch,input=arglist)
}
## if(non.det==TRUE){
## fres <- list(faa=faa[,,-1,drop=F],naa=naa[,,-1,drop=F],biom=biom[,,-1,drop=F],
## ssb=ssb[,,-1,drop=F],wcaa=wcaa[,,-1,drop=F],caa=caa[,,-1,drop=F],
## M=M[,,-1,drop=F],rps=rps.mat[,-1,drop=F],
## maa=maa[,,-1,drop=F],vbiom=apply(biom[,,-1,drop=F],c(2,3),sum,na.rm=T),
## eaa=eaa[,-1,drop=F],
## waa=waa[,,-1,drop=F],waa.catch=waa.catch[,,-1,drop=F],currentF=currentF,
## vssb=apply(ssb[,,-1,drop=F],c(2,3),sum,na.rm=T),vwcaa=vwcaa[,-1,drop=F],
## years=fyears,fyear.year=fyear.year,ABC=ABC,recfunc=recfunc,rec.arg=rec.arg,
## waa.year=waa.year,maa.year=maa.year,multi=multi,multi.year=multi.year,
## Frec=Frec,rec.new=rec.new,pre.catch=pre.catch,input=arglist)
## }
class(fres) <- "future"
invisible(fres)
}
forward.calc.mat2 <- function(fav,nav,Mv,plus.group=TRUE){
nage <- max(which(!is.na(nav[,1])))#length(fav)
na.age <- which(is.na(nav[-1,1]))
# naa <- matrix(NA,nage,dim(nav)[[2]])
naa <- matrix(NA,dim(nav)[[1]],dim(nav)[[2]])
# for(a in 2:(nage-1)){
naa[-c(1,nage,na.age),] <- nav[-c(nage,nage-1,na.age),]*
exp(-fav[-c(nage,nage-1,na.age),]-Mv[-c(nage,nage-1,na.age),])
# }
naa[nage,] <- nav[nage-1,]*exp(-fav[nage-1,]-Mv[nage-1,])
pg <- nav[nage,]*exp(-fav[nage,]-Mv[nage,])
if(plus.group) naa[nage,] <- naa[nage,] + pg
return(naa)
}
caa.est.mat <- function(naa,saa,waa,M,catch.obs,Pope){
saa <- saa/max(saa)
tmpfunc <- function(logx,catch.obs=catch.obs,naa=naa,saa=saa,waa=waa,M=M,out=FALSE,Pope=Pope){
x <- exp(logx)
if(isTRUE(Pope)){
caa <- naa*(1-exp(-saa*x))*exp(-M/2)
}
else{
caa <- naa*(1-exp(-saa*x-M))*saa*x/(saa*x+M)
}
wcaa <- caa*waa
if(out==FALSE){
return(log((sum(wcaa,na.rm=T)-catch.obs)^2))
}
else{
return(caa)
}
}
tmp <- optimize(tmpfunc,log(c(0.000001,10)),catch.obs=catch.obs,naa=naa,saa=saa,waa=waa,M=M,Pope=Pope,out=FALSE)#,tol=.Machine$double.eps)
tmp2 <- tmpfunc(logx=tmp$minimum,catch.obs=catch.obs,naa=naa,saa=saa,waa=waa,M=M,Pope=Pope,out=TRUE)
return(list(x=exp(tmp$minimum),caa=tmp2))
}
# HS用; ARには対応していないが、残差リサンプリングには対応している
HS.rec <- function(ssb,vpares,#deterministic=FALSE,
rec.resample=NULL,
rec.arg=list(a=1000,b=1000,sd=0.1,
resample=FALSE,resid=0, # 残差リサンプリングする場合、resample=TRUEにして、residにリサンプリングする残差(対数)を入れる
bias.correction=TRUE)){
rec0 <- ifelse(ssb>rec.arg$b,rec.arg$a*rec.arg$b,rec.arg$a*ssb)
if(!isTRUE(rec.arg$resample)){
if(isTRUE(rec.arg$bias.correction)){
rec <- rec0*exp(rnorm(length(ssb),-0.5*(rec.arg$sd)^2,rec.arg$sd))
}
else{
rec <- rec0*exp(rnorm(length(ssb),0,rec.arg$sd))
}
}
else{
if(isTRUE(rec.arg$bias.correction)){
rec <- c(rec0[1],
exp(log(rec0[-1])+sample(rec.arg$resid,length(ssb)-1,replace=TRUE))/mean(exp(rec.arg$resid))
)
}
else{
rec <- c(rec0[1],exp(log(rec0[-1])+sample(rec.arg$resid,length(ssb)-1,replace=TRUE)))
}
}
return(list(rec=rec,rec.resample=rec.arg$resid)) # 暫定的変更
}
# RI用; ARには対応していないが、残差リサンプリングには対応している
RI.rec <- function(ssb,vpares,#deterministic=FALSE,
rec.resample=NULL,
rec.arg=list(a=1000,b=1000,sd=0.1,
resample=FALSE,resid=0, # 残差リサンプリングする場合、resample=TRUEにして、residにリサンプリングする残差(対数)を入れる
bias.correction=TRUE)){
rec0 <- rec.arg$a*ssb*exp(-rec.arg$b*ssb) # rec.arg$a*ssb/(1+rec.arg$b*ssb)
# rec0 <- ifelse(ssb>rec.arg$b,rec.arg$a*rec.arg$b,rec.arg$a*ssb)
if(!isTRUE(rec.arg$resample)){
if(isTRUE(rec.arg$bias.correction)){
rec <- rec0*exp(rnorm(length(ssb),-0.5*(rec.arg$sd)^2,rec.arg$sd))
}
else{
rec <- rec0*exp(rnorm(length(ssb),0,rec.arg$sd))
}
}
else{
if(isTRUE(rec.arg$bias.correction)){
rec <- c(rec0[1],exp(log(rec0[-1])+sample(rec.arg$resid,length(ssb)-1,replace=TRUE))/mean(exp(rec.arg$resid)))
}
else{
rec <- c(rec0[1],exp(log(rec0[-1])+sample(rec.arg$resid,length(ssb)-1,replace=TRUE)))
}
}
return(list(rec=rec,rec.resample=rec.arg$resid)) # 暫定的変更
}
# RI用; ARには対応していないが、残差リサンプリングには対応している
BH.rec <- function(ssb,vpares,#deterministic=FALSE,
rec.resample=NULL,
rec.arg=list(a=1000,b=1000,sd=0.1,
resample=FALSE,resid=0, # 残差リサンプリングする場合、resample=TRUEにして、residにリサンプリングする残差(対数)を入れる
bias.correction=TRUE)){
rec0 <- rec.arg$a*ssb/(1+rec.arg$b*ssb)
# rec0 <- ifelse(ssb>rec.arg$b,rec.arg$a*rec.arg$b,rec.arg$a*ssb)
if(!isTRUE(rec.arg$resample)){
if(isTRUE(rec.arg$bias.correction)){
rec <- rec0*exp(rnorm(length(ssb),-0.5*(rec.arg$sd)^2,rec.arg$sd))
}
else{
rec <- rec0*exp(rnorm(length(ssb),0,rec.arg$sd))
}
}
else{
if(isTRUE(rec.arg$bias.correction)){
rec <- c(rec0[1],
exp(log(rec0[-1])+sample(rec.arg$resid,length(ssb)-1,replace=TRUE))/mean(exp(rec.arg$resid))
)
}
else{
rec <- c(rec0[1],exp(log(rec0[-1])+sample(rec.arg$resid,length(ssb)-1,replace=TRUE)))
}
}
return(list(rec=rec,rec.resample=rec.arg$resid)) # 暫定的変更
}
# リサンプリング用(HS.rec, BH.rec, RI.recと同等)
resample.rec <- function(ssb,vpares,#deterministic=FALSE,
rec.resample=NULL,
rec.arg=list(a=1000,b=1000,sd=0.1,
resid=0, # 残差リサンプリングする場合、resample=TRUEにして、residにリサンプリングする残差(対数)を入れる
SR="HS",# or "BH","RI"
bias.correction=TRUE)){
if(rec.arg$SR=="HS") rec0 <- ifelse(ssb>rec.arg$b,rec.arg$a*rec.arg$b,rec.arg$a*ssb)
if(rec.arg$SR=="BH") rec0 <- rec.arg$a*ssb/(1+rec.arg$b*ssb)
if(rec.arg$SR=="RI") rec0 <- rec.arg$a*ssb*exp(-rec.arg$b*ssb)
if(!isTRUE(rec.arg$resample)){
if(isTRUE(rec.arg$bias.correction)){
rec <- rec0*exp(rnorm(length(ssb),-0.5*(rec.arg$sd)^2,rec.arg$sd))
}
else{
rec <- rec0*exp(rnorm(length(ssb),0,rec.arg$sd))
}
}
else{
if(isTRUE(rec.arg$bias.correction)){
rec <- c(rec0[1],exp(log(rec0[-1])+sample(rec.arg$resid,length(ssb)-1,replace=TRUE))/mean(exp(rec.arg$resid)))
}
else{
rec <- c(rec0[1],exp(log(rec0[-1])+sample(rec.arg$resid,length(ssb)-1,replace=TRUE)))
}
}
return(list(rec=rec,rec.resample=rec.arg$resid)) # 暫定的変更
}
# しきい値を設定し、その境界前後でリサンプリングする残差を変える
resample_2block.rec <- function(ssb,vpares,#deterministic=FALSE,
rec.resample=NULL,
rec.arg=list(a=1000,b=1000,sd=0.1,
resid.lower=0, # しきい値よりも小さいときの残差のセット
resid.higher=0, # しきい値よりも大きいときの残差のセット
ssb.threshold=0, # しきい値
SR="HS", # or "BH", "RI"
bias.correction=TRUE)){
if(rec.arg$SR=="HS") rec0 <- ifelse(ssb>rec.arg$b,rec.arg$a*rec.arg$b,rec.arg$a*ssb)
if(rec.arg$SR=="BH") rec0 <- rec.arg$a*ssb/(1+rec.arg$b*ssb)
if(rec.arg$SR=="RI") rec0 <- rec.arg$a*ssb*exp(-rec.arg$b*ssb)
mean.bias.lower <- 1/mean(exp(rec.arg$resid.lower))
mean.bias.higher <- 1/mean(exp(rec.arg$resid.higher))
resid.set <- rep(0,length(ssb[-1]))
tmp <- ssb[-1]<rec.arg$ssb.threshold
if(isTRUE(rec.arg$bias.correction)){
mean.bias <- rep(1,length(ssb[-1]))
mean.bias[tmp] <- mean.bias.lower
mean.bias[!tmp] <- mean.bias.higher
resid.set[tmp] <- sample(rec.arg$resid.lower,sum(tmp),replace=TRUE)
resid.set[!tmp] <- sample(rec.arg$resid.higher,sum(!tmp),replace=TRUE)
det.bias <- ifelse(ssb[1]<rec.arg$ssb.threshold,
mean.bias.lower,mean.bias.higher)
# cat(det.bias,"\n")
# cat(mean(mean.bias),"\n")
rec <- c(rec0[1],#*det.bias,
exp(log(rec0[-1])+resid.set)*mean.bias)
# if(det.bias<1) browser()
}
else{
resid.set[tmp] <- sample(rec.arg$resid.lower,sum(tmp),replace=TRUE)
resid.set[!tmp] <- sample(rec.arg$resid.higher,sum(!tmp),replace=TRUE)
rec <- c(rec0[1],exp(log(rec0[-1])+resid.set))
}
# if(isTRUE(rec.arg$bias.correction)){
# mean.bias <- mean(exp(c(rec.arg$resid.lower,rec.arg$resid.higher)))
# rec <- c(rec0[1],exp(log(rec0[-1])+resid.set))/mean.bias
# }
# else{
# rec <- c(rec0[1],exp(log(rec0[-1])+resid.set))
# }
return(list(rec=rec,rec.resample=rec.arg$resid)) # 暫定的変更
}
# Hockey-stick(bias.correctionのオプションは削除。どうせするので)
HS.recAR <- function(ssb,vpares,#deterministic=FALSE,
rec.resample=NULL,
rec.arg=list(a=1000,b=1000,#gamma=0.01,
sd=0.1, rho=0,
resid=0)#, bias.correction=TRUE)
){
## 再生産関係からの予測値
# rec0 <- rec.arg$a*(ssb+sqrt(rec.arg$b^2+(rec.arg$gamma^2)/4)-sqrt((ssb-rec.arg$b)^2+(rec.arg$gamma^2)/4))
rec0 <- ifelse(ssb>rec.arg$b,rec.arg$a*rec.arg$b,rec.arg$a*ssb)
rec <- rec0*exp(rec.arg$rho*rec.arg$resid) # 自己相関込みの予測値
rec <- rec*exp(rnorm(length(ssb),-0.5*rec.arg$sd2^2,rec.arg$sd))
new.resid <- log(rec/rec0)+0.5*rec.arg$sd2^2
return(list(rec=rec,rec.resample=new.resid))
}
# Beverton-Holt
BH.recAR <- function(ssb,vpares,deterministic=FALSE,rec.resample=NULL,
rec.arg=list(a=1000,b=1000,sd=0.1,bias.correction=TRUE)){
rec0 <- rec.arg$a*ssb/(1+rec.arg$b*ssb)
rec <- rec0*exp(rec.arg$rho*rec.arg$resid) # 自己相関込みの予測値
rec <- rec*exp(rnorm(length(ssb),-0.5*rec.arg$sd2^2,rec.arg$sd))
new.resid <- log(rec/rec0)+0.5*rec.arg$sd2^2
return(list(rec=rec,rec.resample=new.resid))
}
# Ricker
RI.recAR <- function(ssb,vpares,deterministic=FALSE,rec.resample=NULL,
rec.arg=list(a=1000,b=1000,sd=0.1,bias.correction=TRUE)){
rec0 <- rec.arg$a*ssb*exp(-rec.arg$b*ssb)
rec <- rec0*exp(rec.arg$rho*rec.arg$resid) # 自己相関込みの予測値
rec <- rec*exp(rnorm(length(ssb),-0.5*rec.arg$sd2^2,rec.arg$sd))
new.resid <- log(rec/rec0)+0.5*rec.arg$sd2^2
return(list(rec=rec,rec.resample=new.resid))
}
## Allee effect (depensation) ありの再生産関係
# Hockey-stick
HS.recAR2 <- function(ssb,vpares,#deterministic=FALSE,
rec.resample=NULL,
rec.arg=list(a=1000,b=1000,
sd=0.1, rho=0, c=1,
resid=0)#, bias.correction=TRUE)
){
## 再生産関係からの予測値
# rec0 <- rec.arg$a*(ssb+sqrt(rec.arg$b^2+(rec.arg$gamma^2)/4)-sqrt((ssb-rec.arg$b)^2+(rec.arg$gamma^2)/4))
a <- rec.arg$a
b <- rec.arg$b
c <- rec.arg$c
rec0 <- ifelse(ssb>b,b*a,a*b*(ssb/b)^c)
# rec0 <- ifelse(ssb>rec.arg$b,rec.arg$a*rec.arg$b,rec.arg$a*ssb)
rec <- rec0*exp(rec.arg$rho*rec.arg$resid) # 自己相関込みの予測値
rec <- rec*exp(rnorm(length(ssb),-0.5*rec.arg$sd2^2,rec.arg$sd))
new.resid <- log(rec/rec0)+0.5*rec.arg$sd2^2
return(list(rec=rec,rec.resample=new.resid))
}
# Beverton-Holt
BH.recAR2 <- function(ssb,vpares,deterministic=FALSE,rec.resample=NULL,
rec.arg=list(a=1000,b=1000,sd=0.1,rho=0,c=1,bias.correction=TRUE)){
a <- rec.arg$a
b <- rec.arg$b
c <- rec.arg$c
rec0 <- (a/b)/(1+1/(b*ssb)^c)
# rec0 <- rec.arg$a*ssb/(1+rec.arg$b*ssb)
rec <- rec0*exp(rec.arg$rho*rec.arg$resid) # 自己相関込みの予測値
rec <- rec*exp(rnorm(length(ssb),-0.5*rec.arg$sd2^2,rec.arg$sd))
new.resid <- log(rec/rec0)+0.5*rec.arg$sd2^2
return(list(rec=rec,rec.resample=new.resid))
}
# Ricker
RI.recAR2 <- function(ssb,vpares,deterministic=FALSE,rec.resample=NULL,
rec.arg=list(a=1000,b=1000,sd=0.1,rho=0,c=1,bias.correction=TRUE)){
a <- rec.arg$a
b <- rec.arg$b
c <- rec.arg$c
rec0 <- a/(b*exp(1))*(b*ssb)^c*exp(c*(1-b*ssb))
# rec0 <- rec.arg$a*ssb*exp(-rec.arg$b*ssb)
rec <- rec0*exp(rec.arg$rho*rec.arg$resid) # 自己相関込みの予測値
rec <- rec*exp(rnorm(length(ssb),-0.5*rec.arg$sd2^2,rec.arg$sd))
new.resid <- log(rec/rec0)+0.5*rec.arg$sd2^2
return(list(rec=rec,rec.resample=new.resid))
}
plot.futures <- function(fres.list,conf=c(0.1,0.5,0.9),target="SSB",legend.text="",xlim.tmp=NULL,y.scale=1){
if(target=="SSB") aa <- lapply(fres.list,function(x) apply(x$vssb[,-1],1,quantile,probs=conf))
if(target=="Biomass") aa <- lapply(fres.list,function(x) apply(x$vbiom[,-1],1,quantile,probs=conf))
if(target=="Catch") aa <- lapply(fres.list,function(x) apply(x$vwcaa[,-1],1,quantile,probs=conf))
if(target=="Recruit"){
if(is.null(x$recruit)) x$recruit <- x$naa
aa <- lapply(fres.list,function(x) apply(x$recruit[,-1],1,quantile,probs=conf))
}
if(is.null(xlim.tmp)) xlim.tmp <- as.numeric(range(unlist(sapply(aa,function(x) colnames(x)))))
plot(0,max(unlist(aa)),type="n",xlim=xlim.tmp,
ylim=y.scale*c(0,max(unlist(aa))),xlab="Year",ylab=target)
lapply(1:length(aa),function(i) matpoints(colnames(aa[[i]]),t(aa[[i]]),col=i,type="l",lty=c(2,1,2)))
legend("bottomright",col=1:length(aa),legend=legend.text,lty=1)
invisible(aa)
}
plot.future <- function(fres0,ylim.tmp=NULL,xlim.tmp=NULL,vpares=NULL,what=c(TRUE,TRUE,TRUE),conf=0.1,N.line=0,
label=c("Biomass","SSB","Catch"),is.legend=TRUE,add=FALSE,col=NULL,...){
## 暗黙に、vssbなどのmatrixの1列目は決定論的なランの結果と仮定している
if(is.null(col)) col <- 1
matplot2 <- function(x,add=FALSE,...){
if(add==FALSE) matplot(rownames(x),x,type="l",lty=c(2,1,2),col=col,xlab="Year",...)
if(add==TRUE) matpoints(rownames(x),x,type="l",lty=c(2,1,2),col=col,xlab="Year",...)
}
if(is.null(xlim.tmp)) xlim.tmp <- range(as.numeric(rownames(fres0$vssb)))
if(what[1]){
matplot2(x <- t(apply(fres0$vbiom[,-1],1,quantile,probs=c(conf,0.5,1-conf))),
ylim=c(0,ifelse(is.null(ylim.tmp),max(x),ylim.tmp[1])),
xlim=xlim.tmp,
ylab=label[1],main=label[1],add=add,...)
points(rownames(fres0$vbiom),apply(fres0$vbiom[,-1],1,mean),type="b",pch=1)
points(rownames(fres0$vbiom),as.numeric(fres0$vbiom[,1]),type="b",pch=3)
if(!is.null(vpares)){
points(colnames(vpares$baa),colSums(vpares$baa),type="o",pch=20)
}
if(N.line>0) matpoints(rownames(fres0$vbiom),fres0$vbiom[,2:(N.line+1)],col="gray",type="l",lty=1)
}
if(what[2]){
matplot2(x <- t(apply(fres0$vssb[,-1],1,quantile,probs=c(conf,0.5,1-conf))),
ylim=c(0,ifelse(is.null(ylim.tmp),max(x),ylim.tmp[2])),
xlim=xlim.tmp,
ylab=label[2],main=label[2],add=add,...)
points(rownames(fres0$vssb),apply(fres0$vssb[,-1],1,mean),type="b",pch=1)
points(rownames(fres0$vssb),as.numeric(fres0$vssb[,1]),type="b",pch=3)
if(!is.null(fres0$input$Frec))
if(!is.null(fres0$input$Frec$scenario))
if(fres0$input$Frec$scenario!="catch.mean"){
abline(h=fres0$input$Frec$Blimit,col=2)
abline(v=fres0$input$Frec$future.year,col=2)
}
if(!is.null(vpares)){
points(colnames(vpares$ssb),colSums(vpares$ssb),type="o",pch=20)
}
if(N.line>0) matpoints(rownames(fres0$vssb),fres0$vssb[,2:(N.line+1)],col="gray",type="l",lty=1)
}
if(what[3]){
matplot2(x <- t(apply(fres0$vwcaa[,-1],1,quantile,probs=c(conf,0.5,1-conf))),
ylim=c(0,ifelse(is.null(ylim.tmp),max(x),ylim.tmp[3])),
xlim=xlim.tmp,
ylab=label[3],main=label[3],add=add,...)
points(rownames(fres0$vwcaa),apply(fres0$vwcaa[,-1],1,mean),type="b",pch=1)
points(rownames(fres0$vwcaa),as.numeric(fres0$vwcaa[,1]),type="b",pch=3)
if(!is.null(fres0$input$Frec))
if(fres0$input$Frec$scenario=="catch.mean"){
abline(h=fres0$input$Frec$Blimit,col=2)
abline(v=fres0$input$Frec$future.year,col=2)
}
if(!is.null(vpares)){
points(colnames(vpares$baa),colSums(vpares$input$dat$caa*vpares$input$dat$waa),type="o",pch=20)
}
if(N.line>0) matpoints(rownames(fres0$vwcaa),fres0$vwcaa[,2:(N.line+1)],col="gray",type="l",lty=1)
}
if(is.legend){
if(sum(what)>1) plot(1:10,type = "n",ylab = "", xlab = "", axes = F)
legend("topleft",lty=c(NA,NA,1,2),legend=c("Deterministic","Mean","Median",paste(100-(conf*2)*100,"%conf")),pch=c(3,1,NA,NA))
}
}
#print.future <- function(fres){ # S3 method を使いたいんですが、まだいまいちわかりません
# cat(fres$ABC[1])
#}
#
ref.F2 <- function(res0,target.year=c(2018,2023),current.year=2011,Blim,
interval=c(0,3),...){
ssb <- apply(res0$ssb,2,sum)
Frec <- numeric()
Frec[1] <- ssb[current.year]/Blim
for(i in 1:length(target.year)){
tmpfunc <- function(x,res0,Blim,...){
fres <- future.vpa(res0=res0,multi=x,...)
cat(x," ")
return((fres$vssb[rownames(fres$vssb)==target.year[i]]-Blim)^2)
}
Frec[i+1] <- optimize(tmpfunc,interval=interval,res0=res0,Blim=Blim,...)$minimum
}
return(Frec)
}
# 2012. 8. 3 -- 管理基準値計算は外に出す
getABC <- function(res.vpa, # VPAの結果
res.ref, # 管理基準値計算の結果
res.future, # 将来予測計算の結果
ref.case="all",
multi=NULL,
N=NULL,
SSBcur=1000,
Blim=1000,Bban=0,
target.year=NULL, # NULLの場合,ABC.year+4
catch.year=NULL, # 2013:2017など、漁獲量の平均を出したい期間、NULLの場合、ABC.year:ABC.year+4
is.plot=TRUE){
if(all(ref.case=="all")) ref.case <- names(res.ref$summary)
if(all(is.null(multi))) multi <- rep(1,length(ref.case))
nref <- length(ref.case)
ABC.year <- res.future$input$ABC.year
if(is.null(target.year)) target.year <- ABC.year+4
ABC <- wariai <- aveF <- catch5u <- catch5l <- upperSSBlim <- upperSSBcur <- SSBlim <- SSBcur.tmp <- rep(NA,nref)
names(ABC) <- names(wariai) <- names(aveF) <- paste(ref.case,"x",round(multi,3))
wcatch <- matrix(NA,5,nref,dimnames=list(((min(ABC.year)):(min(ABC.year)+4)),names(aveF)))
fres <- list()
i.tmp <- match(ref.case,names(res.ref$summary))
if(any(is.na(i.tmp)))
stop(paste("ref.case specification of is wrong!"))
years <- res.future$year
currentF <- res.ref$Fcurrent["max"] * res.ref$sel
N <- ifelse(is.null(N),dim(res.future$naa)[[3]],N)
for(i in 1:nref){
tmp <- res.ref$summary[i.tmp[i]][1,1] * res.ref$sel
tmp <- max(tmp,na.rm=T)/max(currentF,na.rm=T)*multi[i]
tmpF <- tmp * currentF
aveF[i] <- mean(tmpF,na.rm=T)
input.tmp <- res.future$input
input.tmp$multi <- tmp
input.tmp$is.plot <- FALSE
input.tmp$N <- N
# Frecで使われたシードはとっておかないといけない=> seedはFrecの引数の外に出すこと!
input.tmp$Frec <- NULL
fres[[i]] <- do.call(future.vpa, input.tmp)
ABC[i] <- fres[[i]]$ABC[1]
# browser()
if(res.future$input$ts>1){ # ts>2の場合、漁獲量などの計算は暦年を使う
input.tmp <- res.future$input
input.tmp$multi <- tmp
input.tmp$ts <- 1
input.tmp$is.plot <- FALSE
input.tmp$ABC.year <- ABC.year <- floor(min(input.tmp$ABC.year))
input.tmp$waa <- input.tmp$maa <- input.tmp$M <- input.tmp$waa.catch <- NULL
input.tmp$N <- N
fres[[i]] <- do.call(future.vpa, input.tmp)
years <- fres[[i]]$year
}
wariai[i] <- sum(fres[[i]]$wcaa[,years==ABC.year,1],na.rm=T)/
sum(fres[[i]]$biom[,years==ABC.year,1],na.rm=T)
catch.year <- (ABC.year):(ABC.year+4)
wcatch[,i] <- apply(fres[[i]]$vwcaa[years %in% (catch.year),-1],1,mean,na.rm=T)
catch5u[i] <- quantile(fres[[i]]$vwcaa[years==max(catch.year),-1],probs=0.9) # catchは2017年
catch5l[i] <- quantile(fres[[i]]$vwcaa[years==max(catch.year),-1],probs=0.1)
tmp.year <- years %in% target.year
if(is.null(SSBcur)) SSBcur <- fres[[i]]$vssb[years==(ABC.year),1]
SSBcur.tmp[i] <- SSBcur
upperSSBlim[i] <- sum(fres[[i]]$vssb[tmp.year,-1]>Blim)/N*100 # SSBは2018年当初まで
upperSSBcur[i] <- sum(fres[[i]]$vssb[tmp.year,-1]>SSBcur)/N*100
SSBlim[i] <- Blim
}
if(is.plot){
par(mfrow=c(1,2),mar=c(4,4,2,1))
vssb <- apply(res.vpa$ssb,2,sum,na.rm=T)/1000
x <- sapply(fres,function(x) x$vssb[,1])/1000
plot(range(c(as.numeric(names(vssb)),years)),
c(0,max(x)*1.1),type="n",xlab="Year",ylab="SSB (x1000)")
matpoints(years,x,col=1:nref,type="l",lty=1,
ylim=c(0,max(x)))
points(as.numeric(names(vssb)),vssb,type="b")
abline(h=c(SSBlim/1000,SSBcur/1000),col="gray")
title("SSB in deterministic runs")
plot(0,axes=F,xlab="",ylab="")
legend("topleft",col=1:nref,lty=1,legend=names(ABC))
}
average <- apply(wcatch,2,mean)
res.ref$ABC <- rbind(aveF,wariai,catch5l,catch5u,average,
upperSSBcur,SSBcur.tmp,upperSSBlim,SSBlim,ABC)
rownames(res.ref$ABC)[3] <- paste("catch5l during ",min(catch.year),"-",max(catch.year),sep="")
rownames(res.ref$ABC)[4] <- paste("catch5u during ",min(catch.year),"-",max(catch.year),sep="")
rownames(res.ref$ABC)[5] <- paste("average catch during ",min(catch.year),"-",max(catch.year),sep="")
rownames(res.ref$ABC)[6] <- paste("upperSSBcur at",target.year)
rownames(res.ref$ABC)[8] <- paste("upperSSBlim at",target.year)
fres0 <- fres
write.table(round(res.ref$ABC,2),sep="\t")
save(fres0,file="fres0.R") # 将来予測の全結果はfres0.Rにてセーブされている
# Kobe chartの作成
kobe.array <- array(NA,dim=c(length(fres),nrow(fres[[1]]$vssb),5))
dimnames(kobe.array) <- list(names(ABC),rownames(fres[[1]]$vssb),
c("catch","Biomass","SSB","upperBlimit","upperBban"))
for(i in 1:length(fres)){
kobe.array[i,,] <- as.matrix(get.kobematrix(fres[[i]],
Blim=Blim,Bban=Bban,ssb=TRUE))
}
return(list(ABC=res.ref$ABC,kobe.matrix=kobe.array))
}
#----------------------------------------------------------------------
#---------- 加入に関する関数。魚種specific -------------------
#----------------------------------------------------------------------
#-------------- VPA mode 用関数 -------------------
caa.est <- function(naa,saa,waa,M,catch.obs,Pope){
saa <- saa/max(saa)
tmpfunc <- function(x,catch.obs=catch.obs,naa=naa,saa=saa,waa=waa,M=M,out=FALSE,Pope=Pope){
if(isTRUE(Pope)){
caa <- naa*(1-exp(-saa*x))*exp(-M/2)
}
else{
caa <- naa*(1-exp(-saa*x-M))*saa*x/(saa*x+M)
}
wcaa <- caa*waa
if(out==FALSE){
return((sum(wcaa,na.rm=T)-catch.obs)^2)
}
else{
return(caa)
}
}
tmp <- optimize(tmpfunc,c(0,5),catch.obs=catch.obs,naa=naa,saa=saa,waa=waa,M=M,Pope=Pope,out=FALSE)
tmp2 <- tmpfunc(x=tmp$minimum,catch.obs=catch.obs,naa=naa,saa=saa,waa=waa,M=M,Pope=Pope,out=TRUE)
return(list(x=tmp$minimum,caa=tmp2))
}
#---------------- 結果の確かめ用関数 ---------------------
# --------USAGE-------
# tdata <- get.tdata("vpa_results.csv")
# check.res(res.pms,list(fres,fres),tdata,digits=2,type="%")
get.data <- function(tfile){
tmpdata <- read.csv(tfile,header=F,as.is=F,colClasses="character")
flags <- which(substr(tmpdata[,1],1,1)=="#")
tlist <- list()
for(i in 1:(length(flags)-1)){
tmp <- tmpdata[(flags[i]+1):(flags[i+1]-1),]
if(dim(tmp)[[1]]>1){
dimnames(tmp) <- list(tmp[,1],tmp[1,])
tmp <- tmp[,!apply(tmp=="",2,all)]
tlist[[i]] <- sapply((tmp[-1,-1]),as.numeric)
}
else{
tlist[[i]] <- as.numeric(tmp[tmp!=""])
}
}
names(tlist)[1:4] <- c("naa","faa","Biomass","Fc.at.age")
dimnames(tlist[[3]])[[1]] <- c("SSB","Biomass")
for(i in 1:tlist[[5]]){
names(tlist)[(4+(i-1)*4+1):(4+(i*4))] <- c("fnaa","ffaa","fwcaa","ABC")
}
return(tlist)
}
### 結果の入出力
## 結果の出力
out.vpa <- function(res=NULL, # VPA result
rres=NULL, # reference point
fres=NULL, # future projection result (not nessesarily)
ABC=NULL,
filename="vpa" # filename without extension
){
old.par <- par()
exit.func <- function(){
# par(old.par)
dev.off()
options(warn=0)
}
on.exit(exit.func())
csvname <- paste(filename,".csv",sep="")
pdfname <- paste(filename,".pdf",sep="")
pdf(pdfname)
par(mfrow=c(3,2),mar=c(3,3,2,1))
options(warn=-1)
write.table2 <- function(x,title.tmp="",is.plot=TRUE,...){
if(is.plot){
if(!is.null(dim(x))){
matplot(colnames(x),t(x),type="b",ylim=c(0,max(x,na.rm=T)),pch=substr(rownames(x),1,1))
}
else{
barplot(x)
}
title(title.tmp)
}
if(!is.null(dim(x))){
tmp <- matrix("",nrow(x)+1,ncol(x)+1)
tmp[-1,-1] <- as.character(unlist(x))
tmp[-1,1] <- rownames(x)
tmp[1,-1] <- colnames(x)
}
else{
tmp <- x
}
write.table(tmp,append=T,sep=",",quote=FALSE,file=csvname,col.names=F,row.names=F,...)
}
write(paste("# RVPA outputs at ",date()," & ",getwd()),file=csvname)
if(!is.null(res)){
write("# VPA results",file=csvname, append=T)
write("\n# catch at age",file=csvname,append=T)
write.table2(res$input$dat$caa,title.tmp="Catch at age")
write("\n# maturity at age",file=csvname,append=T)
write.table2(res$input$dat$maa,title.tmp="Maturity at age")
write("\n# weight at age for biomass calculation",file=csvname,append=T)
write.table2(res$input$dat$waa,title.tmp="Weight at age (for biomass)")
if(!is.null(res$input$dat$waa.catch)){
write("\n# weight at age for catch calculation",file=csvname,append=T)
write.table2(res$input$dat$waa.catch,title.tmp="Weight at age (for catch)")
}
write("\n# M at age",file=csvname,append=T)
write.table2(res$input$dat$M,title.tmp="M at age")
write("\n# fishing mortality at age",file=csvname,append=T)
write.table2(res$faa,title.tmp="F at age")
write("\n# Current F",file=csvname,append=T)
write.table2(res$Fc.at.age,title.tmp="Current F")
write("\n# numbers at age",file=csvname,append=T)
write.table2(res$naa,title.tmp="Numbers at age")
write("\n# total and spawning biomass ",file=csvname,append=T)
x <- rbind(colSums(res$ssb),colSums(res$baa),colSums(res$wcaa))
rownames(x) <- c("Spawning biomass","Total biomass","Catch biomass")
write.table2(x,title.tmp="Total and spawning biomass")
}
if(!is.null(rres)){
write("\n# Reference points",file=csvname,append=T)
write.table2(rres$summary,title.tmp="Future F at age",is.plot=F)
}
if(!is.null(fres)){
write("\n# future projection results",file=csvname,append=T)
write("\n# future F at age",file=csvname,append=T)
write.table2(fres$faa[,,1],title.tmp="Future F at age")
write("\n# future numbers at age",file=csvname,append=T)
write.table2(fres$naa[,,1],title.tmp="Future numbers at age")
write("\n# future total and spawning biomass",file=csvname,append=T)
x <- rbind(fres$vssb[,1],fres$vbiom[,1],fres$vwcaa[,1])
rownames(x) <- c("Spawning biomass","Total biomass","Catch biomass")
write.table2(x,title.tmp="Future total, spawning and catch biomass")
}
if(!is.null(ABC)){
write("\n# ABC summary",file=csvname,append=T)
write.table2(ABC$ABC,title.tmp="Future F at age",is.plot=F)
write("\n# Kobe matrix",file=csvname,append=T)
for(i in 1:dim(ABC$kobe.matrix)[[3]]){
write(paste("\n# ",dimnames(ABC$kobe.matrix)[[3]][i]),
file=csvname,append=T)
write.table2(ABC$kobe.matrix[,,i],
title.tmp=dimnames(ABC$kobe.matrix)[[3]][i],is.plot=T)
}
}
}
###
read.vpa <- function(tfile,
caa.label="catch at age",
maa.label="maturity at age",
waa.label="weight at age",
waa.biomass.label="weight at age for biomass calculation",
waa.catch.label="weight at age for catch calculation",
M.label="M at age",
faa.label="fishing mortality at age",
Fc.label="Current F",
naa.label="numbers at age",
Blimit=NULL,
Pope=NULL, # VPA計算時にどっちを使っているか入れる(TRUE or FALSE)。デフォルトはNULLでcaa,faa,naaの関係から自動判別するが、自動判別の結果はcatで出力されるので、それをみて正しく判断されているか確認してください。
fc.year=NULL){
tmpdata <- read.csv(tfile,header=F,as.is=F,colClasses="character")
tmpfunc <- function(tmpdata,label,type=NULL){
flags <- which(substr(tmpdata[,1],1,1)=="#")
flag.name <- tmpdata[flags,1]
flag.name <- gsub("#","",flag.name)
flag.name <- gsub(" ","",flag.name)
get.flag <- which(flag.name==gsub(" ","",label))
{if(length(get.flag)>0){
tdat <- tmpdata[(flags[get.flag]+1):(flags[min(which(flags[get.flag]<flags))]-1),]
if(!is.null(type)){
tdat <- tdat[,!apply(tdat=="",2,all)]
tdat <- as.numeric(tdat)
}
else{
tmp.names <- list(tdat[-1,1],tdat[1,-1])
tdat <- tdat[,!apply(tdat=="",2,all)]
tdat <- tdat[!apply(tdat=="",1,all),]
tdat <- tdat[,!apply(is.na(tdat),2,all)]
tdat <- tdat[!apply(is.na(tdat),1,all),]
tdat <- sapply((tdat[-1,-1]),as.numeric)
tmp.names <- lapply(tmp.names,function(x) x[x!=""])
tmp.names <- lapply(tmp.names,function(x) x[!is.na(x)])
dimnames(tdat) <- tmp.names
tdat <- as.data.frame(tdat)
}
}
else{
tdat <- NULL
}}
tdat
}
dres <- list()
dres$naa <- tmpfunc(tmpdata,naa.label)
dres$faa <- tmpfunc(tmpdata,faa.label)
dres$Fc.at.age <- tmpfunc(tmpdata,Fc.label,type="Fc")
dres$input <- list()
dres$input$dat <- list()
dres$input$dat$maa <- tmpfunc(tmpdata,maa.label)
dres$input$dat$caa <- tmpfunc(tmpdata,caa.label)
dres$input$dat$M <- tmpfunc(tmpdata,M.label)
dres$input$dat$waa <- tmpfunc(tmpdata,waa.label)
if(is.null(dres$input$dat$waa)) dres$input$dat$waa <- tmpfunc(tmpdata,waa.biomass.label)
dres$input$dat$waa.catch <- tmpfunc(tmpdata,waa.catch.label)
if(is.null(dres$input$dat$waa.catch)) dres$input$dat$waa.catch <- dres$input$dat$waa
dres$ssb <- dres$input$dat$waa * dres$input$dat$maa * dres$naa
dres$ssb <- as.data.frame(dres$ssb)
dres$baa <- dres$input$dat$waa * dres$naa
dres$baa <- as.data.frame(dres$baa)
# setting total catch
dres$wcaa <- dres$input$dat$waa.catch * dres$input$dat$caa
dres$wcaa <- as.data.frame(dres$wcaa)
dres$Blimit <- Blimit
## catch at ageの計算時にpopeの近似式を使っているかどうか、通常は外から情報として与えてほしいところだが、与えられない場合、入力されたcaa,faa,naaの関係を見て、Popeで計算されているのかそうでないのかを判断してdres$input$Popeに入れる
if(is.null(Pope)){
caa.pope <- dres$naa*(1-exp(-dres$faa))*exp(-dres$input$dat$M/2)
diff.pope <- mean(unlist(dres$input$dat$caa/caa.pope))
faa <- dres$faa
M <- dres$input$dat$M
caa.bara <- dres$naa*faa/(faa+M)*(1-exp(-faa-M))
diff.bara <- mean(unlist(dres$input$dat$caa/caa.bara))
if(abs(1-mean(diff.bara))>abs(1-mean(diff.pope))){
dres$input$Pope <- TRUE
cat("Pope is TRUE... OK?\n")
}
else{
dres$input$Pope <- FALSE
cat("Pope is FALSE... OK?\n")
}
}
else{
dres$input$Pope <- Pope
}
if(is.null(dres$Fc.at.age) && !is.null(fc.year)) dres$Fc.at.age <- apply(dres$faa[,colnames(dres$faa)%in%fc.year],1,mean)
return(dres)
}
#type="TorF" # true or false
#type="diff" # excel-RVPA
#type="%" # (excel-RVPA)/excel
check.res <- function(res,fres,tdata,digits=3,type="%"){
check.twomats <- function(mat1,mat2,digits=3,type="%"){
if(!is.null(colnames(mat1))){
tmp1 <- mat1[,colnames(mat1)%in%colnames(mat2)]
tmp2 <- mat2[,colnames(mat2)%in%colnames(mat1)]
}
else{
tmp1 <- mat1
tmp2 <- mat2
}
if(type=="TorF"){
tmp <- round(tmp1,digits) == round(tmp2,digits)
}
if(type=="diff"){
tmp <- round(tmp1-tmp2,digits)
}
if(type=="%"){
tmp <- round((tmp1-tmp2)/tmp1*100,digits)
}
return(tmp)
}
naa.res <- check.twomats(tdata$naa,res$naa,digits=digits,type=type)
faa.res <- check.twomats(tdata$faa,res$faa,digits=digits,type=type)
fcaa.res <- check.twomats(tdata$Fc.at.age,res$Fc.at.age,digits=digits,type=type)
tmp.list <- list(naa=naa.res,faa=faa.res,Fc.at.age=fcaa.res)
return(tmp.list)
}
solv.Feq <- function(cvec,nvec,mvec){
Fres <- rep(0,length(cvec))
# cat(nvec," ")
for(i in 1:length(cvec)){
F0 <- cvec[i]/nvec[i]
F1 <- cvec[i]*(F0+mvec[i])/nvec[i]/(1-exp(-F0-mvec[i]))
if(round(cvec[i],6)<round(nvec[i],6)){
while(abs(F0-F1)>0.0001 ){
F0 <- F1
F1 <- cvec[i]*(F0+mvec[i])/nvec[i]/(1-exp(-F0-mvec[i]))
if(F0-F1==-Inf) cat("\n",cvec[i]," ",nvec[i]," \n")
}
Fres[i] <- F1
}
else{
Fres[i] <- 10
cat("Warning: catch exceeded tot_num at: ",i," ",
round(cvec[i],6)," ",round(nvec[i],6),"\n")
}
}
Fres
}
forward.calc.simple <- function(fav,nav,Mv,plus.group=TRUE){
nage <- length(nav)#length(fav)
naa <- rep(NA,nage)
naa[c(-1,-nage)] <- nav[c(-nage,-(nage-1))]*exp(-fav[c(-nage,-(nage-1))]-Mv[c(-nage,-(nage-1))])
naa[nage] <- nav[nage-1]*exp(-fav[nage-1]-Mv[nage-1])
pg <- nav[nage]*exp(-fav[nage]-Mv[nage])
if(plus.group) naa[nage] <- naa[nage] + pg
return(naa)
}
forward.calc.mat <- function(fav,nav,Mv,plus.group=TRUE){
nage <- max(which(!is.na(nav[,1])))#length(fav)
na.age <- which(is.na(nav[,1]))
# naa <- matrix(NA,nage,dim(nav)[[2]])
naa <- matrix(NA,dim(nav)[[1]],dim(nav)[[2]])
# for(a in 2:(nage-1)){
naa[c(-1,-nage,-na.age),] <- nav[c(-nage,-(nage-1),-na.age),]*
exp(-fav[c(-nage,-(nage-1),-na.age),]-Mv[c(-nage,-(nage-1),-na.age),])
# }
naa[nage,] <- nav[nage-1,]*exp(-fav[nage-1,]-Mv[nage-1,])
pg <- nav[nage,]*exp(-fav[nage,]-Mv[nage,])
if(plus.group) naa[nage,] <- naa[nage,] + pg
return(naa)
}
get.kobematrix <- function(fres,Blim=0,Bban=0,ssb=TRUE){
if(isTRUE(ssb)) tmp <- fres$vssb[,-1]
else tmp <- fres$vbiom[,-1]
res <- data.frame(
# 漁獲量
catch.deterministic=fres$vwcaa[,1],
# 資源量
biom.deterministic=fres$vbiom[,1],
# 親魚量
ssb.deterministic=fres$vssb[,1],
# Blim回復確率
probability.upper.Blim=apply(tmp>Blim,1,mean)*100,
# Bban以上確率
probability.upper.Bban=apply(tmp>Bban,1,mean)*100)
return(res)
}
############
# RVPAの結果からMSYを計算する関数
# 主に使うのはSR.est(再生産関係をフィットし、MSYを計算)とSR.plot(フィットした結果をプロット)
############
############
# 使い方
############
if(0){
# マサバ太平洋のデータを読み込み; modelAはvpaの帰り値
modelA <- readRDS("modelA_res.Rdata")
# MSY計算
res1 <- SR.est(modelA,
what.est=c(TRUE,TRUE,TRUE), # HS,BH,RIのどれをフィットするか。
bref.year=2013:2015, # 生物パラメータを用いる期間
years=c(1970:2013), # 観測されたSR関係を用いる期間
er.log=TRUE, # 誤差。TRUEで対数正規誤差
fc.year=2013:2015, # MSY計算のさいに選択率を平均する期間
seed=1 # 乱数の種。この値を変えると乱数が変わるので結果も変わる
)
res1$summary # 推定パラメータ、管理基準値の確認
# 再生産パラメータa,bはエクセルとほぼ一致するはずだが、管理基準値は確率的シミュレーションをもとに計算しているので、エクセルとは必ずしも一致しない。±5%くらいの違いはあるみたい
# 結果のプロット(HSのみ)
res1.pred <- plot.SR(res1,what.plot=c("hs"))
# 結果のプロット(HS,BH,RIを全て)
res1.pred <- plot.SR(res1,what.plot=c("hs","bh","ri"))
allplot(res1) # 要約表・グラフの出力
}
pred.RI <- function(SSB,a,b) a*SSB*exp(-b*SSB)
pred.BH <- function(SSB,a,b) a*SSB/(1+b*SSB)
pred.HS <- function(SSB,a,b,gamma) a*(SSB+sqrt(b^2+gamma^2/4)-sqrt((SSB-b)^2+gamma^2/4))
pred.SL <- function(SSB,a) a*SSB
##
get.stat <- function(fout,eyear=0,hsp=NULL,tmp.year=NULL){
col.target <- ifelse(fout$input$N==0,1,-1)
tmp <- as.numeric(fout$vssb[(nrow(fout$vssb)-eyear):nrow(fout$vssb),col.target])
lhs <- sum(tmp<hsp)/length(tmp)
if(is.null(tmp.year)) tmp.year <- (nrow(fout$vwcaa)-eyear):nrow(fout$vwcaa)
a <- data.frame("catch.mean"=mean(fout$vwcaa[tmp.year,col.target]),
"catch.sd"=sd(fout$vwcaa[tmp.year,col.target]),
"catch.geomean"=geomean(fout$vwcaa[tmp.year,col.target]),
"catch.median"=median(fout$vwcaa[tmp.year,col.target],na.rm=T),
"catch.det"=mean(fout$vwcaa[tmp.year,1],na.rm=T),
"catch.L10"=quantile(fout$vwcaa[tmp.year,col.target],na.rm=T,probs=0.1),
"catch.H10"=quantile(fout$vwcaa[tmp.year,col.target],na.rm=T,probs=0.9),
"ssb.mean"=mean(fout$vssb[tmp.year,col.target]),
"ssb.sd"=sd(fout$vssb[tmp.year,col.target]),
"ssb.geomean"=geomean(fout$vssb[tmp.year,col.target]),
"ssb.median"=median(fout$vssb[tmp.year,col.target],na.rm=T),
"ssb.det"=mean(fout$vssb[tmp.year,1],na.rm=T),
"ssb.L10"=quantile(fout$vssb[tmp.year,col.target],na.rm=T,probs=0.1),
"ssb.H10"=quantile(fout$vssb[tmp.year,col.target],na.rm=T,probs=0.9),
"biom.mean"=mean(fout$vbiom[tmp.year,col.target]),
"biom.sd"=sd(fout$vbiom[tmp.year,col.target]),
"biom.geomean"=geomean(fout$vbiom[tmp.year,col.target]),
"biom.median"=median(fout$vbiom[tmp.year,col.target],na.rm=T),
"biom.det"=mean(fout$vbiom[tmp.year,1],na.rm=T),
"biom.L10"=quantile(fout$vbiom[tmp.year,col.target],na.rm=T,probs=0.1),
"biom.H10"=quantile(fout$vbiom[tmp.year,col.target],na.rm=T,probs=0.9),
"lower.HSpoint"=lhs,
"Fref2Fcurrent"=fout$multi
)
a$U.mean <- a$catch.mean/a$biom.mean
a$U.median <- a$catch.median/a$biom.median
a$U.geomean <- a$catch.geomean/a$biom.geomean
a$U.det <- a$catch.det/a$biom.det
a$catch.CV <- a$catch.sd/a$catch.mean
a$ssb.CV <- a$ssb.sd/a$ssb.mean
a$biom.CV <- a$biom.sd/a$biom.mean
Faa <- as.data.frame(t(fout$multi * fout$input$res0$Fc.at.age))
colnames(Faa) <- paste("F",dimnames(fout$naa)[[1]],sep="")
a <- cbind(a,Faa)
return(a)
}
get.stat2 <- function(fout,unit.waa=1,eyear=2,hsp=NULL,tmp.year=NULL){
col.target <- ifelse(fout$input$N==0,1,-1)
if(is.null(tmp.year)) tmp.year <- (nrow(fout$vwcaa)-eyear):nrow(fout$vwcaa)
nage <- dim(fout$naa)[[1]]
tb <- fout$naa * fout$waa * unit.waa
if(is.null(fout$waa.catch)) fout$waa.catch <- fout$waa
tc <- fout$caa * fout$waa.catch * unit.waa
ssb <- fout$naa * fout$waa *fout$maa * unit.waa
tb.mat <- tc.mat <- ssb.mat <- matrix(0,nage,6)
for(i in 1:nage){
tb.mat[i,1] <- mean(tb[i,tmp.year,col.target])
tb.mat[i,2] <- median(tb[i,tmp.year,col.target])
tb.mat[i,3] <- geomean(tb[i,tmp.year,col.target])
tb.mat[i,4] <- mean(tb[i,tmp.year,1])
tb.mat[i,5:6] <- quantile(tb[i,tmp.year,col.target],probs=c(0.1,0.9),na.rm=T)
tc.mat[i,1] <- mean(tc[i,tmp.year,col.target])
tc.mat[i,2] <- median(tc[i,tmp.year,col.target])
tc.mat[i,3] <- geomean(tc[i,tmp.year,col.target])
tc.mat[i,4] <- mean(tc[i,tmp.year,1])
tc.mat[i,5:6] <- quantile(tc[i,tmp.year,col.target],probs=c(0.1,0.9),na.rm=T)
ssb.mat[i,1] <- mean(ssb[i,tmp.year,col.target])
ssb.mat[i,2] <- median(ssb[i,tmp.year,col.target])
ssb.mat[i,3] <- geomean(ssb[i,tmp.year,col.target])
ssb.mat[i,4] <- mean(ssb[i,tmp.year,1])
ssb.mat[i,5:6] <- quantile(ssb[i,tmp.year,col.target],probs=c(0.1,0.9),na.rm=T)
}
tc.mat <- as.numeric(tc.mat)
tb.mat <- as.numeric(tb.mat)
ssb.mat <- as.numeric(ssb.mat)
# MA; mean, ME; median, GM; geometric mean
names(tc.mat) <- c(paste("TC-MA-A",1:nage,sep=""),paste("TC-ME-A",1:nage,sep=""),
paste("TC-GM-A",1:nage,sep=""),paste("TC-DE-A",1:nage,sep=""),
paste("TC-L10-A",1:nage,sep=""),paste("TC-H10-A",1:nage,sep=""))
names(tb.mat) <- c(paste("TB-MA-A",1:nage,sep=""),paste("TB-ME-A",1:nage,sep=""),
paste("TB-GM-A",1:nage,sep=""),paste("TB-DE-A",1:nage,sep=""),
paste("TB-L10-A",1:nage,sep=""),paste("TB-H10-A",1:nage,sep=""))
names(ssb.mat) <- c(paste("SSB-GA-A",1:nage,sep=""),paste("SSB-ME-A",1:nage,sep=""),
paste("SSB-GM-A",1:nage,sep=""),paste("SSB-DE-A",1:nage,sep=""),
paste("SSB-L10-A",1:nage,sep=""),paste("SSB-H10-A",1:nage,sep=""))
return(as.data.frame(t(c(tb.mat,tc.mat,ssb.mat))))
}
get.stat3 <- function(fout,eyear=0,hsp=NULL,tmp.year=NULL,unit.waa=1){
col.target <- ifelse(fout$input$N==0,1,-1)
tmp <- as.numeric(fout$vssb[(nrow(fout$vssb)-eyear):nrow(fout$vssb),col.target])
lhs <- sum(tmp<hsp)/length(tmp)
if(is.null(tmp.year)) tmp.year <- (nrow(fout$vwcaa)-eyear):nrow(fout$vwcaa)
a <- data.frame("catch.mean"=mean(fout$vwcaa[tmp.year,col.target]),
"catch.sd"=sd(fout$vwcaa[tmp.year,col.target]),
"catch.geomean"=geomean(fout$vwcaa[tmp.year,col.target]),
"catch.median"=median(fout$vwcaa[tmp.year,col.target],na.rm=T),
"catch.det"=mean(fout$vwcaa[tmp.year,1],na.rm=T),
"catch.L10"=quantile(fout$vwcaa[tmp.year,col.target],na.rm=T,probs=0.1),
"catch.H10"=quantile(fout$vwcaa[tmp.year,col.target],na.rm=T,probs=0.9),
"ssb.mean"=mean(fout$vssb[tmp.year,col.target]),
"ssb.sd"=sd(fout$vssb[tmp.year,col.target]),
"ssb.geomean"=geomean(fout$vssb[tmp.year,col.target]),
"ssb.median"=median(fout$vssb[tmp.year,col.target],na.rm=T),
"ssb.det"=mean(fout$vssb[tmp.year,1],na.rm=T),
"ssb.L10"=quantile(fout$vssb[tmp.year,col.target],na.rm=T,probs=0.1),
"ssb.H10"=quantile(fout$vssb[tmp.year,col.target],na.rm=T,probs=0.9),
"biom.mean"=mean(fout$vbiom[tmp.year,col.target]),
"biom.sd"=sd(fout$vbiom[tmp.year,col.target]),
"biom.geomean"=geomean(fout$vbiom[tmp.year,col.target]),
"biom.median"=median(fout$vbiom[tmp.year,col.target],na.rm=T),
"biom.det"=mean(fout$vbiom[tmp.year,1],na.rm=T),
"biom.L10"=quantile(fout$vbiom[tmp.year,col.target],na.rm=T,probs=0.1),
"biom.H10"=quantile(fout$vbiom[tmp.year,col.target],na.rm=T,probs=0.9),
"lower.HSpoint"=lhs,
"Fref2Fcurrent"=fout$multi
)
a$U.mean <- a$catch.mean/a$biom.mean
a$U.median <- a$catch.median/a$biom.median
a$U.geomean <- a$catch.geomean/a$biom.geomean
a$U.det <- a$catch.det/a$biom.det
a$catch.CV <- a$catch.sd/a$catch.mean
a$ssb.CV <- a$ssb.sd/a$ssb.mean
a$biom.CV <- a$biom.sd/a$biom.mean
Faa <- as.data.frame(t(fout$multi * fout$input$res0$Fc.at.age))
colnames(Faa) <- paste("F",dimnames(fout$naa)[[1]],sep="")
res.stat1 <- cbind(a,Faa) # ここまで、get.stat
agename <- dimnames(fout$naa)[[1]]
nage <- dim(fout$naa)[[1]]
tb <- fout$naa * fout$waa * unit.waa
if(is.null(fout$waa.catch)) fout$waa.catch <- fout$waa
tc <- fout$caa * fout$waa.catch * unit.waa
ssb <- fout$naa * fout$waa *fout$maa * unit.waa
tb.mat <- tc.mat <- ssb.mat <- matrix(0,nage,6)
for(i in 1:nage){
tb.mat[i,1] <- mean(tb[i,tmp.year,col.target])
tb.mat[i,2] <- median(tb[i,tmp.year,col.target])
tb.mat[i,3] <- geomean(tb[i,tmp.year,col.target])
tb.mat[i,4] <- mean(tb[i,tmp.year,1])
tb.mat[i,5:6] <- quantile(tb[i,tmp.year,col.target],probs=c(0.1,0.9),na.rm=T)
tc.mat[i,1] <- mean(tc[i,tmp.year,col.target])
tc.mat[i,2] <- median(tc[i,tmp.year,col.target])
tc.mat[i,3] <- geomean(tc[i,tmp.year,col.target])
tc.mat[i,4] <- mean(tc[i,tmp.year,1])
tc.mat[i,5:6] <- quantile(tc[i,tmp.year,col.target],probs=c(0.1,0.9),na.rm=T)
ssb.mat[i,1] <- mean(ssb[i,tmp.year,col.target])
ssb.mat[i,2] <- median(ssb[i,tmp.year,col.target])
ssb.mat[i,3] <- geomean(ssb[i,tmp.year,col.target])
ssb.mat[i,4] <- mean(ssb[i,tmp.year,1])
ssb.mat[i,5:6] <- quantile(ssb[i,tmp.year,col.target],probs=c(0.1,0.9),na.rm=T)
}
tc.mat <- as.numeric(tc.mat)
tb.mat <- as.numeric(tb.mat)
ssb.mat <- as.numeric(ssb.mat)
# mean, ME; median, geomean; geometric mean
names(tc.mat) <- c(paste("TC-mean-A",agename,sep=""),paste("TC-median-A",agename,sep=""),
paste("TC-geomean-A",agename,sep=""),paste("TC-det-A",agename,sep=""),
paste("TC-L10-A",agename,sep=""),paste("TC-H10-A",agename,sep=""))
names(tb.mat) <- c(paste("TB-mean-A",agename,sep=""),paste("TB-median-A",agename,sep=""),
paste("TB-geomean-A",agename,sep=""),paste("TB-det-A",agename,sep=""),
paste("TB-L10-A",agename,sep=""),paste("TB-H10-A",agename,sep=""))
names(ssb.mat) <- c(paste("SSB-GA-A",agename,sep=""),paste("SSB-median-A",agename,sep=""),
paste("SSB-geomean-A",agename,sep=""),paste("SSB-det-A",agename,sep=""),
paste("SSB-L10-A",agename,sep=""),paste("SSB-H10-A",agename,sep=""))
res.stat2 <- as.data.frame(t(c(tb.mat,tc.mat,ssb.mat)))
res.stat <- cbind(res.stat1,res.stat2)
return(res.stat)
}
geomean <- function(x)
{
ifelse(all(x > 0), exp(mean(log(x))), NA)
}
plot.SR <- function(srres,what.plot=c("hs","bh","ri","sl"),xyscale=c(1.3,1.3),xscale=FALSE,is.legend=TRUE,what.sigma=1,FUN="mean",is.MSYline=TRUE,
pick="SSB_MSY"){
col.tmp <- c(rgb(0.3,0.8,0.3,alpha=0.8),rgb(0.8,0.3,0.3,alpha=0.8),rgb(0.3,0.3,0.8,alpha=0.8))
# xscale=TRUEの場合、B0が再生産関係によって異なってくるので、複数の重ね書きはしないこSと!
# tmp <- which(names(srres)==what.plot)
tmp <- which(names(srres)%in%what.plot)
resid <- list()
if(isTRUE(xscale)){
ssb0 <- srres$summary$B0.ssb.mean1[tmp]
xrange <- seq(from=0,to=ssb0,length=100)
}
else{
ssb0 <- 1
xrange <- seq(from=0,to=xyscale[1]*max(srres$dat$SSB,na.rm=T),length=100)
}
plot(x <- srres$dat$SSB/ssb0,y <- srres$dat$R,type="l",pch=20,xlim=range(xrange/ssb0),col="gray",
ylim=c(0,xyscale[2]*max(y,na.rm=T)),xaxs="i",yaxs="i",xlab=ifelse(!xscale,"Spawning biomass (MT)","SB/SB0"),
ylab="Number of recruits",lwd=1)
points(x,y,type="p",pch=20,col=gray(c(seq(from=0.7,to=0,length=length(x)))))
points(rev(x)[1],rev(y)[1],type="p",pch=20,cex=2.5)
for(i in 1:length(what.plot)){
Bmsy <- srres$summary[pick][which(what.plot[i]==rownames(srres$summary)),]
if(what.plot[i]=="hs"){
points(xpred <- xrange/ssb0,
ypred <- pred.HS(SSB=xrange,
a=srres[what.plot[i]][[1]]$a,b=srres[what.plot[i]][[1]]$b,gamma=srres[what.plot[i]][[1]]$gamma),
type="l",lwd=2,col=col.tmp[i],lty=1)
resid[[i]] <- pred.HS(SSB=x,
a=srres[what.plot[i]][[1]]$a,
b=srres[what.plot[i]][[1]]$b,
gamma=srres[what.plot[i]][[1]]$gamma)
resid[[i]] <- log(y)-log(resid[[i]])
Rmsy <- pred.HS(SSB=Bmsy,
a=srres[what.plot[i]][[1]]$a,b=srres[what.plot[i]][[1]]$b,gamma=srres[what.plot[i]][[1]]$gamma)
}
if(what.plot[i]=="bh"|what.plot[i]=="ri"){
if(what.plot[i]=="bh") tmpfunc <- pred.BH
if(what.plot[i]=="ri") tmpfunc <- pred.RI
points(xpred <- xrange/ssb0,
ypred <- tmpfunc(SSB=xrange,
a=srres[what.plot[i]][[1]]$a,b=srres[what.plot[i]][[1]]$b),type="l",lwd=2,col=col.tmp[i],lty=ifelse(what.plot[i]=="bh",2,3))
resid[[i]] <- tmpfunc(SSB=x,
a=srres[what.plot[i]][[1]]$a,
b=srres[what.plot[i]][[1]]$b)
resid[[i]] <- log(y)-log(resid[[i]])
Rmsy <- tmpfunc(SSB=Bmsy,
a=srres[what.plot[i]][[1]]$a,b=srres[what.plot[i]][[1]]$b)
}
if(what.plot[i]=="sl")
points(xrange/ssb0,pred.SL(SSB=xrange,
a=srres[what.plot[i]][[1]]$a),type="l",lwd=1,col=col.tmp[i])
if(is.MSYline){ #abline(v=Bmsy/ssb0,col=i+1,lty=2)
arrows(Bmsy/ssb0,Rmsy*1.2,Bmsy/ssb0,Rmsy,col=col.tmp[i],lty=1,lwd=2,length=.1)
if(Bmsy/ssb0>rev(xrange)[1]){
ymax <- xyscale[2]*max(y,na.rm=T)
arrows(rev(xrange)[2],ifelse(ymax<Rmsy,ymax*0.8,Rmsy*0.8),
rev(xrange)[2],ifelse(ymax<Rmsy,ymax,Rmsy),col=col.tmp[i],lty=1,lwd=2,length=.1)
}
}
}
neg.LL <- sapply(srres[what.plot],function(x) x$res$value)
k <- sapply(srres[what.plot],function(x) length(x$res$par))
n <- length(srres$dat$R)
AICc <- 2*neg.LL+2*k+2*k*(k+1)/(n-k-1)
if(is.legend){
legend("topright",legend=paste(toupper(what.plot[order(AICc)]),round(AICc[order(AICc)],2)),
col=col.tmp[order(AICc)],lwd=1,title="AICc",bg="white",ncol=3)
}
return(list(AICc=AICc,resid=resid,x=xpred,y=ypred))
}
# MSY計算で仮定されている選択率で漁獲したとき何倍になるか? => 計算時間が、、。
# %SPR?MSYを達成した時のFを%SPR換算
plot.Kobe0 <- function(srres,pickB="",what.plot="hs",plot.history=FALSE){
tmp <- which(names(srres)==what.plot)-3
years <- colnames(srres$vpares$ssb)
y <- srres[what.plot][[1]]$Fhist[[1]]$fmulti/srres[what.plot][[1]]$f.msy
x <- as.numeric(colSums(srres$vpares$ssb))/
srres$summary[pickB][[1]][tmp]
x <- x[y>0.001]
years <- years[y>0.001]
y <- y[y>0.001]
plot(x,y,type="n",xlim=c(0,ifelse(max(x)<3,3,max(x,na.rm=T))),
# ylim=c(0,ifelse(max(y)<3,3,max(y))),
ylim=c(0,4),
pch=c(3,rep(1,length(y)-2),20),col=c(1,rep(1,length(y)-2),2),
cex=c(1,rep(1,length(y)-2),2),ylab="F/Fmsy",xlab="B/Bmsy")
polygon(c(-1,1,1,-1),c(-1,-1,1,1),col="khaki1",border=NA)
polygon(c(1,4,4,1),c(-1,-1,1,1),col="olivedrab2",border=NA)
polygon(c(1,4,4,1),c(1,1,6,6),col="khaki1",border=NA)
polygon(c(-1,1,1,-1),c(1,1,6,6),col="indianred1",border=NA)
axis(side=1:2)
points(x,y,type="o",
pch=c(3,rep(1,length(y)-1),20),
col=c(1,rep(1,length(y)-1),1),
cex=c(1,rep(1,length(y)-1),2),ylab="F/Fmsy",xlab="B/Bmsy")
points(rev(x)[1],rev(y)[1],pch=20)
if(isTRUE(plot.history)){
plot(years,y,type="b",ylab="F/Fmsy",ylim=c(0,max(y)))
abline(h=1)
plot(years,x,type="b",xlab="B/Bmsy",ylim=c(0,max(y)))
abline(h=1)
}
invisible(data.frame(years=years,F=y,B=x))
}
plot.Kobe2 <- get.trend <- function(srres,UBdata=NULL,SR="hs",plot.history=FALSE,is.plot=FALSE,pickU="",pickB="",ylab.tmp="U/Umsy",xlab.tmp="SSB/SSBmsy"){
dres <- srres$vpares
tmp <- which(names(srres)==SR)-3
if(is.null(dres$TC.MT)) dres$TC.MT <- as.numeric(colSums(dres$wcaa))
if(is.null(UBdata)){
U <- data.frame(years=as.numeric(colnames(dres$baa)),
U=as.numeric(dres$TC.MT)/as.numeric(colSums(dres$baa,na.rm=T)))
B <- data.frame(years=as.numeric(colnames(dres$ssb)),
B=as.numeric(colSums(dres$ssb)))
Catch <- data.frame(years=as.numeric(colnames(dres$baa)),
C=as.numeric(dres$TC.MT))
UBdata <- merge(U,B)
UBdata <- merge(UBdata,Catch)
# U <- data.frame(years=as.numeric(ts$YEAR),
# U=as.numeric(ts$"TC-MT")/as.numeric(ts$"TB-MT"))
# UBdata$Umsy <- srres$summary$MSY.U.median2[tmp]
# UBdata$Bmsy <- srres$summary$MSY.ssb.median2[tmp]
UBdata$Umsy <- srres$summary[pickU][tmp,]
UBdata$Bmsy <- srres$summary[pickB][tmp,]
UBdata$Uratio <- UBdata$U/UBdata$Umsy
UBdata$Bratio <- UBdata$B/UBdata$Bmsy
}
if(is.plot){
plot(x <- UBdata$Bratio,
y <- UBdata$Uratio,type="n",xlim=c(0,ifelse(max(x)<2,2,max(x,na.rm=T))),
ylim=c(0,ifelse(max(y,na.rm=T)<3,3,max(y,na.rm=T))),
cex=c(1,rep(1,length(y)-2),3),ylab=ylab.tmp,xlab=xlab.tmp)
polygon(c(-1,1,1,-1),c(-1,-1,1,1),col="khaki1",border=NA)
polygon(c(1,6,6,1),c(-1,-1,1,1),col="olivedrab2",border=NA)
polygon(c(1,6,6,1),c(1,1,6,6),col="khaki1",border=NA)
polygon(c(-1,0.5,0.5,-1),c(1,1,6,6),col="indianred1",border=NA)
polygon(c(0.5,1,1,0.5),c(1,1,6,6),col="tan1",border=NA)
polygon(c(-1,0.5,0.5,-1),c(-1,-1,1,1),col="khaki2",border=NA)
polygon(c(0.5,1,1,0.5),c(-1,-1,1,1),col="khaki1",border=NA)
axis(side=1:2)
# points(x,y,type="o",pch=c(3,rep(1,length(y)-2),20),col=c(1,rep(1,length(y)-2),1),cex=c(1,rep(1,length(y)-2),1.5))
points(x,y,type="l",pch=20,col=1,lwd=1)
points(x,y,type="p",pch=20,col=gray(c(seq(from=0.7,to=0,length=length(x)))),cex=1.2)
points(rev(x)[1],rev(y)[1],type="p",pch=20,cex=2.5)
if(isTRUE(plot.history)){
plot(UBdata$years,y,type="b",ylab="F/Fmsy",ylim=c(0,max(y)))
abline(h=1)
plot(UBdata$years,x,type="b",xlab="B/Bmsy",ylim=c(0,max(y)))
abline(h=1)
}}
invisible(UBdata)
}
show.likeprof <- function(res){
x <- tapply(res$hs$surface$obj,list(res$hs$surface$b,res$hs$surface$a),function(x) x)
image(as.numeric(rownames(x)),as.numeric(colnames(x)),log(x/min(x)),col=rev(heat.colors(12)),ylab="a",xlab="b")
contour(as.numeric(rownames(x)),as.numeric(colnames(x)),log(x/min(x)),add=T,nlevels=10,zlim=c(0,0.3))
points(res$hs$b,res$hs$a)
title("Diagnostics")
}
# 単位はcatch at ageの尾数が100万尾、waaがgの場合、重量の単位がちょうどトンになるようになっている。
plot.info <- function(a,xpos=7){
plot(1:(nrow(a)+2),type="n",ylab="",xlab="",axes=F)
units <- ceiling(-1*log10(a[,2]))
units <- units + 2
units <- ifelse(units<0,0,units)
for(i in 1:nrow(a)){
text(1,nrow(a)-i+2,a[i,1],adj=c(0,1),cex=1)
text(xpos,nrow(a)-i+2,format(round(a[i,2],units[i]),big.mark=",",
scientific=F),adj=c(1,1))
}
}
plotfish <- function(image,x,y,size,scale=1,ysize=1){
# image <- readJPEG("../buri.jpg")
xx <- dim(image)[1]/dim(image)[2]
rasterImage(image,
x-size*xinch(1), y-size*yinch(1)*xx*ysize, x+size*xinch(1), y+size*yinch(1)*xx*ysize)
}
menplot <- function(x,y,line.col=1,...){
polygon(c(x,rev(x)),c(y[,1],rev(y[,2])),...)
if(dim(y)[[2]]>2) points(x,y[,3],type="l",lwd=2,col=line.col)
}
menplot2 <- function(xy,probs=c(0.1,0.9),new=FALSE,xlab=NULL,ylab=NULL,...){
xx <- rownames(xy)
yy <- t(apply(xy,1,quantile,probs=c(0.1,0.9)))
if(isTRUE(new)) matplot(xx,yy,type="n",xlab=xlab,ylab=ylab)
menplot(xx,yy,...)
}
plotyield <- function(res00,int.res=NULL,detail.plot=FALSE){
# par(mfrow=c(2,1))
arg.tmp <- res00$farg
arg.tmp$rec.arg$sd <- 0
arg.tmp$N <- 1
# fout.tmp <- do.call(future.vpa2,arg.tmp)
# average
plot(x <- res00$trace$fmulti,y <- res00$trace$catch.mean,type="n",xlim=c(0,max(x)),
xlab="Multiplier to current F",ylab="Catch weight",ylim=c(0,max(res00$trace$catch.det,y)))
menplot(res00$trace$fmulti,cbind(res00$trace$catch.L10,res00$trace$catch.H10),
col=rgb(210/255,94/255,44/255,0.3),border=NA)
## integrate
if(!is.null(int.res)){
points(x,y <- int.res$yield,lty=2,type="o",lwd=1,col="gray")
points(fmax5 <- x[which.max(y)],y[which.max(y)],pch=20,col="gray")
}
points(x <- res00$trace$fmulti,y <- res00$trace$catch.mean,type="l",xlim=c(0,max(x)),
xlab="Multiplier to current F",ylab="Catch weight",ylim=c(0,max(res00$trace$catch.det,y)))
points(fmax1 <- x[which.max(y)],y[which.max(y)],pch=20,col=1)
if(isTRUE(detail.plot)){
# geomean
points(x <- res00$trace$fmulti,y <- res00$trace$catch.geomean,col=2,type="l",xlim=c(0,2))
points(fmax2 <- x[which.max(y)],y[which.max(y)],pch=20,col=2)
# median
points(x <- res00$trace$fmulti,y <- res00$trace$catch.median,col=3,type="l",xlim=c(0,2))
points(fmax3 <- x[which.max(y)],y[which.max(y)],pch=20,col=3)
}
# deteministic
points(x <- res00$trace$fmulti,y <- res00$trace$catch.det,col=4,
type="l",xlim=c(0,2))
points(fmax4 <- x[which.max(y)],y[which.max(y)],pch=20,col=4)
title("Yield vs. F")
## plot CV of yield
par(new=T)
y <- res00$trace$catch.CV
plot(x,y,type="l",lwd=3,
col=rgb(0.8,0.8,0,0.6),axes=F,xlab="",ylab="",
ylim=c(0,ifelse(max(y,na.rm=T)>1.5,1.5,max(y,na.rm=T))))
axis(side=4)
mtext(side=4,"CV of catch",line=2.5,col=rgb(0.8,0.8,0,0.6),cex=0.8)
### plot SSB
plot(x <- res00$trace$fmulti,y <- res00$trace$ssb.mean,type="n",xlim=c(0,max(x)),
xlab="Relative F (to current F)",ylab="SSB")
menplot(res00$trace$fmulti,cbind(res00$trace$ssb.L10,res00$trace$ssb.H10),
col=rgb(40/255,96/255,163/255,0.3),border=NA)
## integrate
if(!is.null(int.res)){
points(x,y <- int.res$ssb,lty=2,type="o",lwd=1,col="gray")
points(fmax5,y[x==fmax5],pch=20,col="gray")
}
points(x <- res00$trace$fmulti,y <- res00$trace$ssb.mean,type="l",xlim=c(0,max(x)),
xlab="Relative F (to current F)",ylab="SSB")
points(fmax1,y[x==fmax1],pch=20,col=1)
if(isTRUE(detail.plot)){
points(x <- res00$trace$fmulti,y <- res00$trace$ssb.geomean,col=2,type="l",xlim=c(0,2))
points(fmax2,y[x==fmax2],pch=20,col=2)
points(x <- res00$trace$fmulti,y <- res00$trace$ssb.median,col=3,type="l",xlim=c(0,2))
points(fmax3,y[x==fmax3],pch=20,col=3)
}
points(x <- res00$trace$fmulti,y <- res00$trace$ssb.det,
col=4,type="l")
points(fmax4,y[x==fmax4],pch=20,col=4)
title("SSB vs. F")
if(!is.null(int.res)){
legend("topright",lty=c(1,1,1,1,2,NA),col=c(1:4,"gray",NA),legend=c("Simple mean","Geometric mean","Median","Deterministic","Integrate","fill: 80% conf"),bty="n")
}
else{
if(isTRUE(detail.plot)){
legend("topright",lty=c(1,1,1,1,NA),col=c(1:4,NA),
legend=c("Simple mean","Geometric mean","Median","Deterministic","fill: 80% conf"))
}
else{
legend("topright",lty=c(1,1,NA),col=c(c(1,4),NA),
legend=c("Simple mean","Deterministic","fill: 80% conf"))
}
}
#### CV plot
par(new=T)
y <- res00$trace$ssb.CV
plot(x,y,type="l",lwd=3,
col=rgb(0.8,0.8,0,0.6),axes=F,xlab="",ylab="",
ylim=c(0,ifelse(max(y,na.rm=T)>1.5,1.5,max(y,na.rm=T))))
axis(side=4)
mtext(side=4,"CV of SSB",line=2.5,col=rgb(0.8,0.8,0,0.6),cex=0.8)
# points(fout.tmp$multi,fout.tmp$vssb[100,1],pch=4)
}
get.SPR <- function(dres,target.SPR=30,Fmax=10,max.age=Inf){
# Fの歴史的な%SPRを見てみる
# 毎年異なるFや生物パラメータに対して、YPR,SPR、SPR0がどのくらい変わっているのか見る(Rコード例2)
# target.SPRが与えられると、target.SPR(%)として与えた数字に対応するSPR値に対するFの乗数も出力する(与えない場合は30%とする)
dres$ysdata <- matrix(0,ncol(dres$faa),5)
dimnames(dres$ysdata) <- list(colnames(dres$faa),c("perSPR","YPR","SPR","SPR0","F/Ftarget"))
for(i in 1:ncol(dres$faa)){
dres$Fc.at.age <- dres$faa[,i] # Fc.at.ageに対象年のFAAを入れる
if(all(dres$Fc.at.age>0)){
byear <- colnames(dres$faa)[i] # 何年の生物パラメータを使うか
# RVPAのref.F関数でYPRなどを計算。
# 配布している1.3から1.4にアップデートしているので、新しいほうの関数を使うこと(返り値がちょっと違う)
a <- ref.F(dres,waa.year=byear,maa.year=byear,M.year=byear,rps.year=2000:2011,
pSPR=round(target.SPR),
F.range=c(seq(from=0,to=ceiling(max(dres$Fc.at.age,na.rm=T)*Fmax),
length=101),max(dres$Fc.at.age,na.rm=T)),plot=FALSE,max.age=max.age)
# YPRと%SPR
dres$ysdata[i,1:2] <- (as.numeric(rev(a$ypr.spr[which(a$ypr.spr$Frange2Fcurrent==1)[1],2:3])))
# SPR
dres$ysdata[i,3] <- a$spr0*dres$ysdata[i,1]/100
# SPR0
dres$ysdata[i,4] <- a$spr0
# relative F
dres$ysdata[i,5] <- 1/a$summary[3,grep("SPR",colnames(a$summary))][1]
}
else{
break;
}
}
dres$ysdata <- as.data.frame(dres$ysdata)
dres$target.SPR <- target.SPR
return(dres)
}
get.SRdata <- function(vpares,R.dat=NULL,SSB.dat=NULL,years=as.numeric(colnames(vpares$naa))){
# R.datとSSB.datだけが与えられた場合、それを使ってシンプルにフィットする
if(!is.null(R.dat) & !is.null(SSB.dat)){
dat <- data.frame(R=R.dat,SSB=SSB.dat,year=1:length(R.dat))
}
else{
# データの整形
n <- ncol(vpares$naa)
L <- as.numeric(rownames(vpares$naa)[1])
dat <- list()
dat$R <- as.numeric(vpares$naa[1,])
dat$SSB <- as.numeric(colSums(vpares$ssb,na.rm=TRUE))
dat$year <- as.numeric(colnames(vpares$ssb))
# 加入年齢分だけずらす
dat$R <- dat$R[(L+1):n]
dat$SSB <- dat$SSB[1:(n-L)]
dat$year <- dat$year[(L+1):n]
# データの抽出
dat <- as.data.frame(dat)
dat <- dat[dat$year%in%years,]
}
class(dat) <- "SRdata"
return(dat[c("year","SSB","R")])
}
plot.SRdata <- function(SRdata){
plot(SRdata$SSB,SRdata$R,xlab="SSB",ylab="R",xlim=c(0,max(SRdata$SSB)),ylim=c(0,max(SRdata$R)))
}
est.MSY <- function(vpares,farg,
seed=farg$seed,
nyear=NULL,
eyear=0, # 将来予測の最後のeyear+1年分を平衡状態とする
# FUN=median, # 漁獲量の何を最大化するか?
FUN=mean, # 漁獲量の何を最大化するか?
N=1000, # stochastic計算するときの繰り返し回数
onlylower.pgy=FALSE,# PGY計算するとき下限のみ計算する(計算時間省略のため)
optim.method="optimize",
max.target="catch.mean", # method="optimize"以外を使うとき、どの指標を最大化するか。他のオプションとしては"catch.median" (漁獲量のmedianの最大化)
calc.yieldcurve=TRUE, # yield curveを正確に計算するかどうか。TRUEだと計算時間が余計にかかる。FALSEだと、yield curveは正確ではない
Blimit=0,
trace.multi=c(seq(from=0,to=0.9,by=0.1),1,seq(from=1.1,to=2,by=0.1),3:5,7,20,100), # Fmsyを探索したり、Yield curveを書くときにグリッドサーチをするときのFの刻み。Fcurrentに対する乗数。Fが異常に大きい場合、親魚=0になって加入=NA
is.plot=TRUE,
PGY=NULL, # PGY管理基準値を計算するかどうか。計算しない場合はNULLを、計算する場合はc(0.8,0.9,0.95)のように割合を入れる
B0percent=NULL, # B0_XX%の管理基準値を計算するかどうか
Bempirical=NULL, # 特定の親魚量をターゲットにする場合
long.term=20, # 世代時間の何倍年後の状態を平衡状態と仮定するか
GT=NULL, # 世代時間を外から与える場合(世代時間の計算は将来予測で使われる年齢別成熟率・自然死亡係数を使っているが、別のパラメータを与えたい場合など、外で計算してここに入れる)
mY=5, # 自己相関を考慮して管理基準値を計算する場合、平衡状態から何年進めるか
resid.year=0, # ARありの場合、最近年何年分の残差を平均するか
current.resid=NULL # 残差の値を直接入れる場合。上の年数が設定されていてもこちらが設定されたらこの値を使う
){
require(tidyverse)
farg$seed <- seed
### 内部で使うための関数定義
## 最小化のための関数
## シミュレーション回数ぶんの漁獲量のFUN(mean, geomean, median)を最大化するFを選ぶ
msy.objfun <- function(x,f.arg,FUN=FUN,eyear=eyear){
f.arg$multi <- x
fout <- do.call(future.vpa,f.arg)
return(-FUN(fout$vwcaa[(nrow(fout$vwcaa)-eyear):nrow(fout$vwcaa),-1]))
}
trace.func <- function(farg,eyear,hsp=0,trace.N=farg$N,
fmulti=c(seq(from=0,to=0.9,by=0.1),1,seq(from=1.1,to=2,by=0.1),3:5,7,20,100)){
trace.res <- NULL
# ssb.array <- array(0,dim=c(farg$nyear,farg$N+1,length(fmulti)))
farg$outtype <- "FULL"
farg$N <- trace.N
for(i in 1:length(fmulti)){
farg$multi <- fmulti[i]
tmp <- do.call(future.vpa,farg)
# ssb.array[,,i] <- tmp$vssb
tmp2 <- get.stat3(tmp,eyear=eyear,hsp=hsp)
trace.res <- rbind(trace.res,tmp2)
if(tmp2$"ssb.mean"<trace.res$"ssb.mean"[1]/1000){
fmulti <- fmulti[1:i]
break()
}
}
trace.res <- as.data.frame(trace.res)
trace.res$fmulti <- fmulti
return(list(table=trace.res))
}
which.min2 <- function(x){
max(which(min(x)==x))
}
target.func <- function(fout,faa0=NULL,mY=5,N=2,seed=1,eyear=4,p=1,beta=NULL,delta=NULL,Blim=0,Bban=0,sd0=NULL,current.resid=NULL){
farg <- fout$input
last.year <- dim(fout$naa)[[2]]
lag <- as.numeric(dimnames(fout$naa)[[1]])[1]
# if(lag==0) SSB.m <- NULL else SSB.m <- fout$ssb[,last.year-lag,]
SSB.m <- fout$ssb[,last.year-lag,]
ssb0 <- SSB.m
farg$seed <- seed
farg$N <- N
farg$nyear <- mY
farg$naa0 <- p*fout$naa[,last.year,]
farg$eaa0 <- fout$eaa[last.year,]+current.resid
farg$ssb0 <- p*ssb0
farg$faa0 <- faa0
farg$beta <- beta
farg$delta <- delta
farg$Blim <- Blim
farg$Bban <- Bban
farg$start.year <- max(as.numeric(colnames(farg$res0$naa)))+1
farg$ABC.year <- farg$start.year
if(!is.null(sd0)) farg$rec.arg$sd <- sd0
farg$Frec <- NULL
fout <- do.call(future.vpa,farg)
out <- get.stat3(fout,eyear=0,hsp=Blimit)
# out <- cbind(out,get.stat2(fout,eyear=0,hsp=Blimit))
return(list(out,fout))
}
### 関数定義おわり
## 世代時間を計算
if(is.null(GT)){
GT <- Generation.Time(vpares,maa.year=farg$maa.year,
M.year=farg$M.year) # Generation Time
}
if(is.null(nyear)){
nyear <- round(GT*long.term)
}
trace.N <- N
years <- sort(as.numeric(rev(names(vpares$naa))[1:5]))
nY <- nyear+1 # これ必要??
## 引数の調整
b0 <- numeric() # B0
fout <- fout0 <- trace <- Fhist <- fout.HS.5par <- list()
farg.org <- farg.tmp <- farg
farg.tmp$outtype <- "FULL"
farg.tmp$nyear <- nyear
farg.tmp$N <- N
farg.tmp$silent <- TRUE
farg.tmp$is.plot <- FALSE
farg.tmp$ABC.year <- max(years)+1
farg.tmp$add.year <- 1
farg.tmp$det.run <- FALSE
if(!is.null(farg.tmp$pre.catch)){
farg.tmp$pre.catch <- NULL # pre.catchオプションがあるとうまくいかないのでなかったことにする
cat("notice: option \"pre.catch\" is turned off in estimating MSY.\n")
}
if(!is.null(farg.tmp$rec.new)){
farg.tmp$rec.new <- NULL # rec.newプションがあるとうまくいかないのでなかったことにする
cat("notice: option \"rec.new\" is turned off in estimating MSY.\n")
}
# B0の計算
farg.tmp$multi <- 0
fout0 <- do.call(future.vpa,farg.tmp)
B0 <- get.stat3(fout0,eyear=eyear,hsp=Blimit)
# B0 <- cbind(B0,get.stat2(fout0,eyear=eyear,hsp=Blimit))
rownames(B0) <- "B0"
trace <- trace.func(farg.tmp,eyear,hsp=Blimit,fmulti=trace.multi,trace.N=trace.N)
xx <- which.max(trace$table$catch.mean)+c(-1,1)
range.tmp <- trace$table$fmulti[xx]
if(xx[1]==0) range.tmp <- c(0,range.tmp)
if(is.na(range.tmp[2])) range.tmp[2] <- max(trace$table$fmulti)*10
farg.tmp$multi <- 1
cat("Estimating MSY\n")
if(optim.method=="optimize"){
tmp <- optimize(msy.objfun,range.tmp,f.arg=farg.tmp,eyear=eyear,FUN=FUN)
# 壁にあたっている限り続ける
while(sum(round(tmp$minimum,3)==range.tmp)>0){
tmp0 <- round(tmp$minimum,3)==range.tmp
range.tmp <- sort(c(range.tmp[tmp0],
range.tmp[tmp0] -2*(mean(range.tmp) - range.tmp[tmp0])))
range.tmp <- ifelse(range.tmp<0,0,range.tmp)
tmp <- optimize(msy.objfun,range.tmp,f.arg=farg.tmp,eyear=eyear,FUN=FUN)
}
farg.msy <- farg.tmp
farg.msy$multi <- tmp$minimum # Fc.at.a * multiがFmsy
cat("F multiplier=",tmp$minimum,"\n")
fout.msy <- do.call(future.vpa,farg.msy)
fout.msy$input$multi <- fout.msy$multi
if(calc.yieldcurve){
trace$table <- rbind(trace$table,trace.func(farg.msy,eyear,hsp=Blimit,trace.N=trace.N,
fmulti=tmp$minimum+c(-0.025,-0.05,-0.075,0,0.025,0.05,0.075))$table)
trace$table <- trace$table[order(trace$table$fmulti),]
}
F.msy <- fout.msy$input$multi*vpares$Fc.at.age
}
# optimizeでなくgridでやる場合
else{
Fmulti <- seq(from=min(range.tmp),to=max(range.tmp),by=0.01)
trace.tmp <- trace.func(farg.tmp,eyear,hsp=Blimit,fmulti=Fmulti,trace.N=trace.N)
farg.msy <- farg.tmp
farg.msy$multi <- trace.tmp$table$fmulti[which.max(unlist(trace.tmp$table[max.target]))]
cat("F multiplier=",farg.msy$multi,"\n")
fout.msy <- do.call(future.vpa,farg.msy)
trace$table <- rbind(trace$table,trace.tmp$table)
trace$table <- trace$table[order(trace$table$fmulti),]
}
MSY <- get.stat3(fout.msy,eyear=eyear)
# MSY <- cbind(MSY,get.stat2(fout.msy,eyear=eyear))
rownames(MSY) <- "MSY"
# cat(" SSB=",MSY$"ssb.mean","\n")
## PGYの計算
fout.PGY <- list()
PGYstat <- NULL
if(!is.null(PGY)){
s <- 1
for(j in 1:length(PGY)){
cat("Estimating PGY ",PGY[j]*100,"%\n")
ttmp <- trace$table$catch.mean-PGY[j]*MSY$catch.mean
ttmp <- which(diff(sign(ttmp))!=0)
frange.list <- list(trace$table$fmulti[ttmp[1]+0:1],
trace$table$fmulti[ttmp[2]+0:1])
if(isTRUE(onlylower.pgy)) i.tmp <- 2 else i.tmp <- 1:2
for(i in i.tmp){
farg.pgy <- farg.tmp
if(sum(is.na(frange.list[[i]]))>0) frange.list[[i]] <- c(0,300)
farg.pgy$Frec <- list(stochastic=TRUE,
future.year=rev(rownames(fout0$vssb))[1],
Blimit=PGY[j]*MSY$catch.mean,
scenario="catch.mean",Frange=frange.list[[i]])
fout.PGY[[s]] <- do.call(future.vpa,farg.pgy)
fout.PGY[[s]]$input$multi <- fout.PGY[[s]]$multi
PGYstat <- rbind(PGYstat,get.stat3(fout.PGY[[s]]))
if(calc.yieldcurve){
trace$table <- rbind(trace$table,trace.func(farg.msy,eyear,hsp=Blimit,trace.N=trace.N,
fmulti=fout.PGY[[s]]$multi+c(-0.025,-0.05,-0.075,0,0.025,0.05,0.075))$table)
trace$table <- trace$table[order(trace$table$fmulti),]
}
s <- s+1
}
}
# PGYstat <- as.data.frame(t(sapply(fout.PGY,get.stat3,eyear=eyear,hsp=Blimit)))
# PGYstat <- cbind(PGYstat,as.data.frame(t(sapply(fout.PGY,get.stat2,eyear=eyear,hsp=Blimit))))
rownames(PGYstat) <- names(fout.PGY) <- paste("PGY",rep(PGY,each=length(i.tmp)),
rep(c("upper","lower")[i.tmp],length(PGY)),sep="_")
}
else{
PGYstat <- NULL
}
###
## B0_%の計算
fout.B0percent <- list()
B0stat <- NULL
if(!is.null(B0percent)){
for(j in 1:length(B0percent)){
cat("Estimating B0 ",B0percent[j]*100,"%\n")
ttmp <- trace$table$ssb.mean-B0percent[j]*B0$ssb.mean
ttmp <- which(diff(sign(ttmp))!=0)
frange.list <- trace$table$fmulti[ttmp[1]+0:1]
farg.b0 <- farg.tmp
farg.b0$Frec <- list(stochastic=TRUE,
future.year=rev(rownames(fout0$vssb))[1],
Blimit=B0percent[j]*B0$ssb.mean,
scenario="ssb.mean",Frange=frange.list)
fout.B0percent[[j]] <- do.call(future.vpa,farg.b0)
fout.B0percent[[j]]$input$multi <- fout.B0percent[[j]]$multi
B0stat <- rbind(B0stat,get.stat3(fout.B0percent[[j]]))
if(calc.yieldcurve){
trace$table <- rbind(trace$table,trace.func(farg.msy,eyear,hsp=Blimit,trace.N=trace.N,
fmulti=fout.B0percent[[j]]$multi+c(-0.025,-0.05,-0.075,0,0.025,0.05,0.075))$table)
trace$table <- trace$table[order(trace$table$fmulti),]
}
}
rownames(B0stat) <- names(fout.B0percent) <- paste("B0-",B0percent*100,"%",sep="")
}
else{
B0stat <- NULL
}
###
## 特定のSSBを目指す場合
fout.Bempirical <- list()
Bempirical.stat <- NULL
if(!is.null(Bempirical)){
for(j in 1:length(Bempirical)){
cat("Estimating B empirical ",Bempirical[j],"\n")
ttmp <- trace$table$ssb.mean-Bempirical[j]
ttmp <- which(diff(sign(ttmp))!=0)
frange.list <- trace$table$fmulti[ttmp[1]+0:1]
farg.ben <- farg.tmp
farg.ben$Frec <- list(stochastic=TRUE,
future.year=rev(rownames(fout0$vssb))[1],
Blimit=Bempirical[j],
scenario="ssb.mean",Frange=frange.list)
fout.Bempirical[[j]] <- do.call(future.vpa,farg.ben)
fout.Bempirical[[j]]$input$multi <- fout.Bempirical[[j]]$multi
Bempirical.stat <- rbind(Bempirical.stat,get.stat3(fout.Bempirical[[j]]))
if(calc.yieldcurve){
trace$table <- rbind(trace$table,trace.func(farg.msy,eyear,hsp=Blimit,trace.N=trace.N,
fmulti=fout.Bempirical[[j]]$multi+c(-0.025,-0.05,-0.075,0,0.025,0.05,0.075))$table)
trace$table <- trace$table[order(trace$table$fmulti),]
}
}
rownames(Bempirical.stat) <- names(fout.Bempirical) <- paste("Ben-",round(Bempirical),"",sep="")
}
else{
Bempirical.stat <- NULL
}
###
refvalue <- bind_rows(MSY,B0,PGYstat,B0stat,Bempirical.stat) %>% as_tibble %>%
mutate(RP_name=c("MSY","B0",rownames(PGYstat),rownames(B0stat),rownames(Bempirical.stat)),
AR=FALSE)
refvalue <- refvalue %>%
mutate(SSB2SSB0=refvalue$ssb.mean/refvalue$ssb.mean[2])
sumvalue <- refvalue %>% select(RP_name,AR,ssb.mean,SSB2SSB0,biom.mean,U.mean,catch.mean,catch.CV,Fref2Fcurrent)
colnames(sumvalue) <- c("RP_name","AR","SSB","SSB2SSB0","B","U","Catch","Catch.CV","Fref/Fcur")
sumvalue <- bind_cols(sumvalue,refvalue[,substr(colnames(refvalue),1,1)=="F"])
### ARありの場合の管理基準値の計算(平衡状態から5年分進めたときの値)
if(resid.year > 0 && is.null(current.resid)){
current.resid <- mean(rev(fout.msy$input$rec.arg$resid)[1:resid.year])
cat("Residuals of ",resid.year," years are averaged as, ",current.resid,"\n")
}
else{
if(resid.year==0){
current.resid <- 0
}
}
lag <- as.numeric(rownames(fout.msy$naa))[1]
eyear <- mY+(lag > 0)*(lag-1)
MSY2 <- target.func(fout.msy,mY=mY,seed=seed,N=N,eyear=mY,current.resid=current.resid)
B02 <- target.func(fout0,mY=mY,seed=seed,N=N,eyear=mY,current.resid=current.resid)
if(!is.null(PGY)){
PGYstat2 <- lapply(1:length(fout.PGY),
function(x) target.func(fout.PGY[[x]],mY=mY,seed=seed,N=N,eyear=mY,current.resid=current.resid))
}
else{
PGYstat2 <- NULL
}
if(!is.null(B0percent)){
B0stat2 <- lapply(1:length(fout.B0percent),
function(x) target.func(fout.B0percent[[x]],mY=mY,seed=seed,N=N,eyear=mY,current.resid=current.resid)
)
}
else{
B0stat2 <- NULL
}
if(!is.null(Bempirical)){
Bempirical.stat2 <- lapply(1:length(fout.Bempirical),
function(x) target.func(fout.Bempirical[[x]],mY=mY,seed=seed,N=N,eyear=mY,current.resid=current.resid)
)
}
else{
Bempirical.stat2 <- NULL
}
refvalue2 <- bind_rows(MSY2[[1]],B02[[1]],
purrr::map_dfr(PGYstat2,function(x) x[[1]]),
purrr::map_dfr(B0stat2,function(x) x[[1]]),
purrr::map_dfr(Bempirical.stat2,function(x) x[[1]])) %>% as_tibble() %>%
mutate(RP_name=refvalue$RP_name,AR=TRUE)
refvalue2 <- refvalue2 %>%
mutate(SSB2SSB0=refvalue$ssb.mean/refvalue$ssb.mean[2])
sumvalue2 <- refvalue2 %>% select(RP_name,AR,ssb.mean,SSB2SSB0,biom.mean,U.mean,catch.mean,catch.CV,Fref2Fcurrent)
colnames(sumvalue2) <- c("RP_name","AR","SSB","SSB2SSB0","B","U","Catch","Catch.CV","Fref/Fcur")
sumvalue2 <- bind_cols(sumvalue2,refvalue2[,substr(colnames(refvalue2),1,1)=="F"])
ssb.ar.mean <- cbind(apply(MSY2[[2]]$vssb,1,mean),
apply(B02[[2]]$vssb,1,mean),
sapply(PGYstat2,function(x) apply(x[[2]]$vssb,1,mean)),
sapply(B0stat2,function(x) apply(x[[2]]$vssb,1,mean)),
sapply(Bempirical.stat2,function(x) apply(x[[2]]$vssb,1,mean)))
ssb.ar.mean <- sweep(matrix(as.numeric(ssb.ar.mean),nrow(ssb.ar.mean),ncol(ssb.ar.mean)),
2,unlist(sumvalue$SSB),FUN="/")
colnames(ssb.ar.mean) <- rownames(sumvalue$SSB)
### 結果のプロットなど
trace$table <- subset(trace$table,fmulti>0)
if(isTRUE(is.plot)){
# plot of yield curve
par(mfrow=c(1,3),mar=c(4,4,2,1))
plot(trace$table$fmulti,trace$table$"ssb.mean"*1.2,type="n",xlab="Fref/Fcurrent",ylab="SSB")
abline(v=sumvalue$Fref2Fcurrent,col="gray")
text(sumvalue$Fref2Fcurrent,max(trace$table$"ssb.mean")*seq(from=1.1,to=0.8,length=nrow(sumvalue)),rownames(sumvalue))
menplot(trace$table$fmulti,cbind(0,trace$table$"ssb.mean"),col="skyblue",line.col="darkblue")
title("Equiribrium SSB")
plot(trace$table$fmulti,trace$table$"catch.mean",type="n",xlab="Fref/Fcurrent",ylab="Catch")
abline(v=sumvalue$Fref2Fcurrent,col="gray")
menplot(trace$table$fmulti,cbind(0,trace$table$"catch.mean"),col="lightgreen",line.col="darkgreen")
title("Equiribrium Catch (Yield curve)")
# plot of the effect of AR
matplot(ssb.ar.mean,type="b",ylab="SSB_MSY_AR/SSB_MSY",xlab="Years from Equiribrium")
legend("topright",col=1:ncol(ssb.ar.mean),legend=rownames(sumvalue),lty=1:ncol(ssb.ar.mean))
title("plot of the effect of AR")
}
## kobe II matrix
#kobe2 <- array(0,dim=c(dim(trace$array)[[1]],dim(trace$array)[[2]],length(sumvalue$SSb)))
#for(i in 1:length(sumvalue$SSB)){
#tmp <- trace$array > sumvalue$SSB[i]
#kobe2[,,i] <- cbind(kobe2,apply(tmp,c(1,2),mean))
# }
#dimnames(kobe2) <- list()
input.list <- list(B0=fout0$input,
msy=fout.msy$input,
pgy=lapply(fout.PGY,function(x) x$input),
B0percent=lapply(fout.B0percent,function(x) x$input))
allsum <- bind_rows(sumvalue,sumvalue2)
allsum$RP.definition <- NA
allsum$RP.definition[allsum$AR==FALSE&allsum$RP_name=="MSY"] <- "Btarget0"
allsum$RP.definition[allsum$AR==FALSE&allsum$RP_name=="PGY_0.9_lower"] <- "Blow0"
allsum$RP.definition[allsum$AR==FALSE&allsum$RP_name=="PGY_0.6_lower"] <- "Blimit0"
allsum$RP.definition[allsum$AR==FALSE&allsum$RP_name=="PGY_0.1_lower"] <- "Bban0"
output <- list(summary=as.data.frame(as.matrix(sumvalue)),
summaryAR=as.data.frame(as.matrix(sumvalue2)),
summary_tb=allsum,
F.msy=F.msy,
all.stat=as.data.frame(as.matrix(refvalue)),
all.statAR=as.data.frame(as.matrix(refvalue2)),
all.stat_tb=bind_rows(refvalue,refvalue2),
trace=trace$table,input.list=input.list,
ssb.ar.mean=ssb.ar.mean)
output$SPR.msy <- calc_MSY_spr(output)
invisible(output)
}
#### function definition
get.perform <- function(fout0,Blimit=0,longyear=50,smallcatch=0.5,N=NULL,
shortyear=c(3,5,10),tmp.year=NULL){
stat1 <- get.stat(fout0,eyear=0,hsp=Blimit,tmp.year=tmp.year)[c("catch.mean","catch.CV","biom.mean","biom.CV","ssb.mean","lower.HSpoint")]
stat2 <- get.stat2(fout0,eyear=0,tmp.year=tmp.year)
stat2 <- data.frame(t(as.data.frame(strsplit(colnames(stat2),"-"))),value=as.numeric(stat2))
rownames(stat2) <- NULL
# waaによる加重平均年齢&組成
xx <- subset(stat2,X1=="TB" & X2=="MA")
nage <- sum(!is.na(xx$value))
tmp <- c(rep(2,ceiling(nage/3)),rep(3,ceiling(nage/3)))
tmp <- c(rep(1,nage-length(tmp)),tmp)
if(sum(tmp==1)==0 & sum(tmp==2)>1) tmp[1] <- 1
xx$bvalue <- xx$value * fout0$waa[,1,1]
xx$waa <- fout0$waa[,1,1]
large.portion1 <- tapply(xx$bvalue[!is.na(xx$bvalue)],tmp,sum,na.rm=T)
stat1$largefish.nature <- large.portion1[names(large.portion1)==3]/sum(large.portion1)
aage.biom <- sum(xx$bvalue * 0:(length(xx$bvalue)-1))/sum(xx$bvalue)
xx <- subset(stat2,X1=="TC" & X2=="MA")
xx$bvalue <- xx$value * fout0$waa[,1,1]
aage.catch <- sum(xx$bvalue * 0:(length(xx$bvalue)-1))/sum(xx$bvalue)
large.portion2 <- tapply(xx$bvalue[!is.na(xx$bvalue)],tmp,sum,na.rm=T)
stat1$largefish.catch <- large.portion2[names(large.portion2)==3]/sum(large.portion2)
# 漁獲量<0.5平均漁獲量の頻度
if(is.null(tmp.year)) tmp.year <- nrow(fout0$vwcaa)
stat1$catch.safe <- 1/mean(fout0$vwcaa[tmp.year,]<smallcatch*mean(fout0$vwcaa[tmp.year,]))
stat1$catch.safe <- ifelse(stat1$catch.safe>longyear,longyear,stat1$catch.safe)
# 親魚量<Blimitの頻度 → 確率の逆数
stat1$ssb.safe <- 1/stat1$"lower.HSpoint"
stat1$ssb.safe <- ifelse(stat1$ssb.safe>longyear,longyear,stat1$ssb.safe)
# ABC.yearから5年目までの平均累積漁獲量
short.catch <- numeric()
for(i in 1:length(shortyear)){
years <- fout0$input$ABC.year:(fout0$input$ABC.year+shortyear[i])
short.catch[i] <- mean(apply(fout0$vwcaa[rownames(fout0$vwcaa)%in%years,-1],2,sum))
}
names(short.catch) <- paste("short.catch",shortyear,sep="")
short.catch <- as.data.frame(t(short.catch))
# 平衡状態になった年
years <- names(fout0$vssb[,1])[-1]
heikou.diff <- which(diff(fout0$vssb[,1])/fout0$vssb[-1,1]<0.01)
if(length(heikou.diff)>0) stat1$eq.year <- years[min(heikou.diff)] else stat1$eq.year <- Inf
dat <- data.frame(stat1,short.catch,aage.biom=aage.biom,aage.catch=aage.catch,effort=fout0$multi,
waa=as.data.frame(t(fout0$waa[,1,1])),meigara=as.data.frame(t(tmp)))
return(dat)
}
plotRadial <- function(index,base=1,col.tmp=NULL,lwd=2,...){
old.par <- par()
layout(matrix(c(1,2),2,1),heights=c(2,1))
index2 <- sweep(matrix(unlist(index),nrow(index),ncol(index)),2,as.numeric(unlist(index[base,])),FUN="/")
if(is.null(col.tmp)) col.tmp <- brewer.pal(nrow(index2-1),"Set1")
radial.plot(index2,rp.type="p",lwd=lwd,show.grid.labels=FALSE,
labels=colnames(index),
radial.lim=c(0,1.5),clockwise=TRUE,start=1,
line.col=c(NA,col.tmp),
poly.col=c(rgb(40/255,96/255,163/255,0.2),rep(NA,nrow(index2)-1)), # MSYだけ色で塗る
...
)
refname <- rownames(index)
par(mar=c(1,0,1,0))
plot(0,10,type="n",axes=FALSE,ylab="")
legend("topleft",legend=refname,
col=c(rgb(40/255,96/255,163/255,0.2),col.tmp),
ncol=2,lwd=c(10,rep(lwd,length(refname)-1)))
layout(matrix(c(1),1,1),heights=c(1))
par(old.par)
invisible(index2)
}
## 管理基準値を取り出す関数
get.Bref <- function(res,SRfunc="hs",B0=c(0.3),SPR0=c(0.3),HS=c(1,1.3),PGY=c("PGY_0.9_upper_hs","PGY_0.9_lower_hs")){
sumref <- res$summary[rownames(res$summary)==SRfunc,]
refpoints <- list()
## MSY管理基準値をピックアップ
refpoints$BMSY <- sumref$"SSB_MSY"
## B0基準の管理基準値はB0×%
## B0の値はmout$summary$"B0(SSB)"にある。1番目がHSの結果
refpoints$B0per <- sumref$"B0(SSB)"[1] * B0 # B0_10,20,30,35,40%の値
names(refpoints$B0per) <- paste(B0*100,"%",sep="")
## B_HS関連の管理基準値
refpoints$BHS <- sumref$b[1] * HS
names(refpoints$BHS) <- paste("B_HSx",HS,sep="")
## B_PGY関連の管理基準値(HSをもとにしたものはPGY.biom.hsにあります)
x <- res$PGY.biom.hs["ssb.mean"]
refpoints$BPGY <- x[match(PGY,rownames(x)),1]
names(refpoints$BPGY) <- PGY
## SSB_current
refpoints$SSBcur <- rev(as.numeric(res$vpares$ssb))[1]
## SSB_max
refpoints$SSBmax <- max(as.numeric(res$vpares$ssb))
return(unlist(refpoints))
}
plot.RP <- function(rdata,RP=NULL,biomass.scale=1,ymax=1,is.text=TRUE){
n <- length(rdata)
rdata <- sort(rdata)
if(is.null(RP)) RP <- names(rdata)
ymax <- ymax * seq(from=0.5,to=1,length=n)
for(j in 1:n){
abline(v=rdata[j]/biomass.scale,lty=1,lwd=2,col=rgb(40/255,96/255,40/255,0.5))
if(isTRUE(is.text)){
text(rdata[j]/biomass.scale,ymax[j],
paste(RP[j],"=\n",format(round(rdata[j]/biomass.scale),big.mark=","),"",sep=""),adj=0)
}
}
}
#### 資源量の上積みグラフを書く
plotBfish <- function(res0, # SR.estの結果
Bref,
unit.waa=1,ssb.max=Inf,
target="hs",biomass.scale=1000,pngfile="fish.png"){
summary <- res0$summary[rownames(res0$summary)==target,]
res00 <- res0[names(res0)==target][[1]]
tres0 <- res00$trace[[1]]
ssb <- res00$trace[[1]]$ssb.mean/biomass.scale
tmp <- substr(colnames(tres0),1,5)=="TB-MA"
tb <- tres0[,tmp]/biomass.scale * unit.waa
tb2 <- sapply(1:ncol(tb),function(x) apply(tb[,1:x,drop=F],1,sum,na.rm=T))
tmp <- substr(colnames(tres0),1,5)=="TC-MA"
tc <- tres0[,tmp]/biomass.scale * unit.waa
tc2 <- sapply(1:ncol(tc),function(x) apply(tc[,1:x,drop=F],1,sum,na.rm=T))
library(png)
if(file.exists(pngfile)) image <- readPNG(pngfile)
else image <- NULL
year.tmp <- rev(colnames(res0$vpares$ssb))[1:5]
range1 <- range(res0$vpares$ssb)/biomass.scale
range2 <- range(as.data.frame(res0$vpares$ssb)[as.character(year.tmp)])/biomass.scale
col.tmp1 <- rgb(40/255,96/255,163/255,seq(from=0.1,to=0.9,length=ncol(tc)))
col.tmp2 <- rgb(210/255,94/255,44/255,seq(from=0.1,to=0.9,length=ncol(tc)))
### plot of SSB
ssb.max <- min(ssb.max,
max(c(range1,summary$"SSB_MSY"),na.rm=T)) *1.5 /biomass.scale
tb3 <- tb2[which(ssb<ssb.max),]
matplot(ssb,tb2,type="n",ylab=paste("Biomass (",biomass.scale," MT)",sep=""),xaxs="i",yaxs="i",
xlab="SSB",
ylim=c(0,max(tb2[which(ssb<ssb.max),])*1.2),xlim=c(0,ssb.max))
# ylim=c(0,max(tb2)),xlim=c(0,ssb.max))
# menplot(range1,cbind(c(-100,-100),rep(max(tb2)*1.5,2)),col=gray(0.9),border=NA)
# menplot(range2,cbind(c(-100,-100),rep(max(tb2)*1.5,2)),col=gray(0.7),border=NA)
# 管理基準値のプロット
plot.RP(Bref,biomass.scale=biomass.scale,ymax=max(tb3)*1.1)
# 過去の時系列
# matpoints(ssb,tb2[,1],type="l",lwd=2,col="gray",lty=3)
points(x <- colSums(res0$vpares$ssb)/biomass.scale,
y <- colSums(res0$vpares$baa)/biomass.scale,type="o",
col=gray(c(seq(from=0.7,to=0,length=length(x)))),pch=20,cex=1.2,
lwd=3)
text(x[1],y[1],colnames(x)[1],adj=0)
text(rev(x[1]),rev(y)[1],rev(colnames(x))[1],adj=0)
## 積み上げグラフ
non.na <- !is.na(ssb)
for(i in 1:ncol(tb2)) menplot(ssb[non.na], cbind(0,tb2)[non.na,i:(i+1)],col=col.tmp1[i],border=NA)
title("Total biomass",line=-1,adj=0.1)
# browser()
## abline(v=summary$"SSB_MSY"/biomass.scale,lty=2)
## abline(v=summary$"Blimit"/biomass.scale,lty=2)
## abline(v=summary$"SSB_HS"/biomass.scale,lty=2)
## text(summary$"SSB_MSY"/biomass.scale,max(tb3)*1.1,
## paste("SSB_MSY=",format(round(summary$"SSB_MSY"/biomass.scale),big.mark=","),"",sep=""),adj=0)
## text(summary$"Blimit"/biomass.scale,max(tb3)*1.0,
## paste("SSB_limit=",format(round(summary$"Blimit"/biomass.scale),big.mark=","),
## "",sep=""),adj=0)
## text(summary$"SSB_HS"/biomass.scale,max(tb3)*1.05,
## paste("SSB_HS=",format(round(summary$"SSB_HS"/biomass.scale),big.mark=","),
## "",sep=""),adj=0)
## catch
if(!is.null(res0$vpares$wcaa)) wcatch <- as.numeric(colSums(res0$vpares$wcaa))
else{
wcatch <- as.numeric(colSums(res0$vpares$input$dat$caa * res0$vpares$input$dat$waa,na.rm=T))*unit.waa
}
matplot(ssb,tc2,type="n",,xaxs="i",yaxs="i",ylab=paste("Catch (",biomass.scale,") MT",sep=""),
xlab="SSB",
# ylim=c(0,max(tc2,wcatch)*1.2),xlim=c(0,ssb.max))
ylim=c(0,max(tc2)*1.2),xlim=c(0,ssb.max))
# menplot(range1,cbind(c(-100,-100),rep(max(tb2)*1.5,2)),col=gray(0.9),border=NA)
# menplot(range2,cbind(c(-100,-100),rep(max(tb2)*1.5,2)),col=gray(0.7),border=NA)
points(x <- as.numeric(colSums(res0$vpares$ssb))/biomass.scale,
y <- wcatch/biomass.scale,pch=20,lwd=3,
type="o",col=gray(c(seq(from=0.7,to=0,length=length(x)))))
plot.RP(Bref,biomass.scale=biomass.scale,ymax=max(tc2)*1.1,is.text=FALSE)
# scale <- max(ssb)/max(tc2) * 0.8
for(i in 1:ncol(tc2)) menplot(ssb[non.na], cbind(0,tc2)[non.na,i:(i+1)],col=col.tmp2[i],border=NA)
## abline(v=summary$"SSB_MSY"/biomass.scale,lty=2)
## abline(v=summary$"Blimit"/biomass.scale,lty=2)
## abline(v=summary$"SSB_HS"/biomass.scale,lty=2)
## text(x[1],y[1],colnames(res0$vpares$ssb)[1],adj=0)
## text(rev(x)[1],rev(y)[1],rev(colnames(res0$vpares$ssb))[1],adj=0)
# points(x <- apply(res00$fout[[1]]$vssb,1,mean)[1:10]/biomass.scale,
# y <- apply(res00$fout[[1]]$vwcaa,1,mean)[1:10]/biomass.scale,col=2,
# type="o",pch=20,lwd=3)
# text(rev(x)[1],rev(y)[1],
# paste("Projection ",rownames(res00$fout[[1]]$vssb)[10],"(F_MSY)",sep=""),adj=-0.1,col=2)
# points(x <- apply(fout0$vssb,1,mean)[1:10]/biomass.scale,
# y <- apply(fout0$vwcaa,1,mean)[1:10]/biomass.scale,col="blue",type="o",pch=20,lwd=3)
# text(rev(x)[1],rev(y)[1],paste("現在のFでの10年将来予測"),adj=-0.1,col="blue")
# 魚のプロット
waa.tmp <- (apply(res0$vpares$input$dat$waa,1,mean))^(1/3)*10
waa.tmp <- waa.tmp/max(waa.tmp) * 0.9
x <- tc2[which.min(abs(ssb-ssb.max*0.88)),]
if(!is.null(image)){
plotfish(image,x=rep(ssb.max*0.88,ncol(tc2)),y=x-diff(c(0,x))/2,
size=waa.tmp*0.8,scale=scale,ysize=1)
}
text(rep(ssb.max*0.9,ncol(tc2)),x-diff(c(0,x))/2,
paste(0:(ncol(tc2)-1),"y/o: ",round(apply(res0$vpares$input$dat$waa,1,mean),0)," g"),cex=1)
title("Total catch",line=-1,adj=0.1)
## 努力量やCVのプロット
tmp <- round(ssb*biomass.scale)>0 & !is.na(ssb)
matplot(ssb,tres0$fmulti,type="l",ylab="Efforts (Current=1)",col=1,lwd=2,
xaxs="i",yaxs="i",xlab=paste("SSB (",biomass.scale,"MT)",sep=""),xlim=c(0,ssb.max),
ylim=c(0,max(tres0$fmulti[tmp]*1.2)))
# menplot(range1,cbind(c(-100,-100),rep(max(tb2)*1.5,2)),col=gray(0.9),border=NA)
# menplot(range2,cbind(c(-100,-100),rep(max(tb2)*1.5,2)),col=gray(0.7),border=NA)
# menplot(ssb[tmp],cbind(0,tres0$fmulti[tmp]),col=rgb(221/255,159/255,33/255,0.5),border=NA)
plot.RP(Bref,biomass.scale=biomass.scale,ymax=max(tc2)*1.1,is.text=FALSE)
title("Efforts",line=-1.5,font=2,adj=0.1)
par(new=T)
# y <- res00$trace[[1]]$ssb.CV
y <- res00$trace[[1]]$catch.CV
plot(ssb,y,type="l",lwd=2,col=2,axes=F,xlab="",ylab="",
ylim=c(0,ifelse(max(y,na.rm=T)>1.5,1.5,max(y,na.rm=T))))
# points(ssb,y,type="l",lwd=2,col=3)
axis(side=4)
mtext(side=4,"Catch CV",line=3,col=2,cex=0.8)
}
## kobe.matrixの計算
# Pr(B<Btarget)のみを返す単純なやつ
get.kobemat <- function(fout,N=fout$input$N,nyear=fout$input$nyear,Btarget=0,
fmulti=seq(from=0.3,to=1,by=0.1)){
multi.org <- 1
fres.short <- list()
farg <- fout$input
farg$Frec <- NULL
farg$N <- N
farg$nyear <- nyear
for(i in 1:length(fmulti)){
farg$multi <- multi.org * fmulti[i]
fres.short[[i]] <- do.call(future.vpa,farg)
}
prob.btarget <- sapply(fres.short,function(x) apply(x$vssb>Btarget,1,mean))
colnames(prob.btarget) <- fmulti
invisible(prob.btarget)
}
# Btargetはベクトルで入力、平均親魚量なども出力
get.kobemat2 <- function(fout,N=fout$input$N,nyear=fout$input$nyear,Btarget=0,
fmulti=seq(from=0.3,to=1,by=0.1),target.name=1:length(Btarget)){
multi.org <- 1
fres.short <- list()
farg <- fout$input
farg$Frec <- NULL
farg$N <- N
farg$nyear <- nyear
for(i in 1:length(fmulti)){
farg$multi <- multi.org * fmulti[i]
fres.short[[i]] <- do.call(future.vpa,farg)
}
# 結果の取り出し
prob.btarget <- list()
for(i in 1:length(Btarget)){
prob.btarget[[i]] <- sapply(fres.short,function(x) apply(x$vssb>Btarget,1,mean))
colnames(prob.btarget[[i]]) <- fmulti
}
names(prob.btarget) <- target.name
# SSB, biomass, catch
ssb <- sapply(fres.short,function(x) apply(x$vssb,1,mean))
biom <- sapply(fres.short,function(x) apply(x$vbiom,1,mean))
catch <- sapply(fres.short,function(x) apply(x$vwcaa, 1,mean))
colnames(ssb) <- colnames(biom) <- colnames(catch) <- fmulti
invisible(list(prob.btarget=prob.btarget,ssb=ssb,biom=biom,catch=catch))
}
plot.kobemat <- function(xx,title.name="",line=0){
yy <- as.data.frame.table(xx)
yy$pch <- 20
yy$color <- "gray"
yy$color[as.numeric(yy[,3])>0.5] <- "red"
yy$color[yy[,3]<0.5 & yy[,3]>0.4] <- "pink"
plot(as.numeric(as.character(yy[,1])),
as.numeric(as.character(yy[,2])),type="n",col=yy$color,xlab="Years",pch=yy$pch,cex=3,
ylab="multiplier to F_current")
abline(h=seq(from=0,to=10,by=0.1),v=2000:2100,col="gray")
points(as.numeric(as.character(yy[,1])),
as.numeric(as.character(yy[,2])),col=yy$color,pch=yy$pch,cex=3)
title(title.name)
abline(h=line,col="red")
text(2+min(as.numeric(as.character(yy[,1]))),line,paste("F_",title.name,sep=""))
legend("topleft",col=c("gray","pink","red"),legend=c("<40%","40-50%",">50%"),
title="Pr(B>Btarget)",pch=20,pt.cex=3,bg="white")
}
plot.kobemat2 <- function(yy,...){
xx <- yy$prob.btarget
for(i in 1:length(xx)){
plot.kobemat(xx[[i]],title.name=names(xx)[i],line=-1)
}
matplot(yy$ssb,type="l",ylim=c(0,max(yy$ssb)),lty=1)
title("SSB",line=-1)
matplot(yy$biom,type="l",ylim=c(0,max(yy$biom)),lty=1)
title("Biomass",line=-1)
matplot(yy$catch,type="l",ylim=c(0,max(yy$catch)),lty=1)
title("Catch",line=-1)
legend("bottomright",col=1:ncol(yy$ssb),legend=colnames(yy$ssb),lty=1,title="Fcurrentx")
}
####################
### 西嶋加筆 # 2018/06/07
## 加入の残差の自己相関を考慮した再生産関係の推定
## L1ノルム(最小絶対値)も推定できる (sigmaはSD)
## TMB = TRUEでmarginal likelihood (.cppファイルが必要)
fit.SR <- function(SRdata,
SR="HS",
method="L2",
AR=1,TMB=FALSE,
hessian=FALSE,w=rep(1,length(SRdata$year)),
length=20,
max.ssb.pred=1.3, # 予測値を計算するSSBの最大値(観測された最大値への乗数)
p0=NULL,
out.AR = FALSE #自己相関係数rhoを外で推定するか
){
argname <- ls()
arglist <- lapply(argname,function(xx) eval(parse(text=xx)))
names(arglist) <- argname
if (AR==0) out.AR <- FALSE
rec <- SRdata$R
ssb <- SRdata$SSB
N <- length(rec)
NN <- sum(w) #likelihoodを計算するサンプル数
# if (SR=="HS") SRF <- function(x,a,b) a*(x+sqrt(b^2+gamma^2/4)-sqrt((x-b)^2+gamma^2/4))
if (SR=="HS") SRF <- function(x,a,b) ifelse(x>b,b*a,x*a)
if (SR=="BH") SRF <- function(x,a,b) a*x/(1+b*x)
if (SR=="RI") SRF <- function(x,a,b) a*x*exp(-b*x)
obj.f <- function(a,b,rho){
resid <- sapply(1:N,function(i) log(rec[i]) - log(SRF(ssb[i],a,b)))
resid2 <- NULL
for (i in 1:N) {
resid2[i] <- ifelse(i==1, resid[i], resid[i]-rho*resid[i-1])
}
if (method == "L2") {
rss <- w[1]*resid2[1]^2*(1-rho^2)
for(i in 2:N) rss <- rss + w[i]*resid2[i]^2
sd <- sqrt(rss/NN)
sd2 <- c(sd/sqrt(1-rho^2), rep(sd,N-1))
obj <- -sum(w*dnorm(resid2,0,sd2,log=TRUE))
} else {
rss <- w[1]*abs(resid2[1])*sqrt(1-rho^2)
for(i in 2:N) rss <- rss + w[i]*abs(resid2[i])
sd <- sum(abs(w*resid2))/NN
sd2 <- c(sd/sqrt(1-rho^2), rep(sd,N-1))
obj <- -sum(w*sapply(1:N, function(i){-log(2*sd2[i])-abs(resid2[i]/sd2[i])}))
}
return(obj)
}
if (is.null(p0)) {
a.range <- range(rec/ssb)
b.range <- range(1/ssb)
if (SR == "HS") b.range <- range(ssb)
grids <- as.matrix(expand.grid(
seq(a.range[1],a.range[2],len=length),
seq(b.range[1],b.range[2],len=length)
))
init <- as.numeric(grids[which.min(sapply(1:nrow(grids),function(i) obj.f(grids[i,1],grids[i,2],0))),])
init[1] <- log(init[1])
init[2] <- ifelse (SR == "HS",-log(max(0.000001,(max(ssb)-min(ssb))/max(init[2]-min(ssb),0.000001)-1)),log(init[2]))
if (AR != 0 && !isTRUE(out.AR)) init[3] <- 0
} else {
init = p0
}
if (SR == "HS") {
if (AR == 0 || out.AR) {
obj.f2 <- function(x) obj.f(exp(x[1]),min(ssb)+(max(ssb)-min(ssb))/(1+exp(-x[2])),0)
} else {
obj.f2 <- function(x) obj.f(exp(x[1]),min(ssb)+(max(ssb)-min(ssb))/(1+exp(-x[2])),1/(1+exp(-x[3])))
}
} else {
if (AR == 0 || out.AR) {
obj.f2 <- function(x) obj.f(exp(x[1]),exp(x[2]),0)
} else {
obj.f2 <- function(x) obj.f(exp(x[1]),exp(x[2]),1/(1+exp(-x[3])))
}
}
opt <- optim(init,obj.f2)
opt <- optim(opt$par,obj.f2,method="BFGS",hessian=hessian)
Res <- list()
Res$input <- arglist
Res$opt <- opt
a <- exp(opt$par[1])
b <- ifelse(SR=="HS",min(ssb)+(max(ssb)-min(ssb))/(1+exp(-opt$par[2])),exp(opt$par[2]))
rho <- ifelse(AR==0,0,ifelse(out.AR,0,1/(1+exp(-opt$par[3]))))
resid <- sapply(1:N,function(i) log(rec[i]) - log(SRF(ssb[i],a,b)))
resid2 <- NULL
for (i in 1:N) {
resid2[i] <- ifelse(i == 1,resid[i], resid[i]-rho*resid[i-1])
}
if (method=="L2") {
rss <- w[1]*resid2[1]^2*(1-rho^2)
for(i in 2:N) rss <- rss + w[i]*resid2[i]^2
sd <- sqrt(rss/NN)
} else {
rss <- w[1]*abs(resid2[1])*sqrt(1-rho^2)
for(i in 2:N) rss <- rss + w[i]*abs(resid2[i])
sd <- sum(abs(w*resid2))/NN
sd <- sqrt(2)*sd
}
# sd <- ifelse(method=="L2",sqrt(sum(w*resid2^2)/(NN-rho^2)),sqrt(2)*sum(abs(w*resid2))/(NN-rho^2))
Res$resid <- resid
Res$resid2 <- resid2
Res$pars <- c(a,b,sd,rho)
if (method!="L2") {
if (AR!=0) {
arres <- ar(resid,aic=FALSE,order.max=1)
Res$pars[3] <- ifelse(arres$ar<0,sd,sqrt(arres$var.pred))
Res$pars[4] <- ifelse(arres$ar<0,0,arres$ar)
}
}
if (AR==1 && out.AR) {
arres <- ar(resid,aic=FALSE,order.max=1,demean=FALSE,method="mle")
Res$pars[3] <- sqrt(arres$var.pred)
Res$pars[4] <- as.numeric(arres$ar)
Res$resid2[2:length(Res$resid2)] <- arres$resid[-1]
Res$AIC.ar <- ar(resid,order.max=1,demean=FALSE,method="mle")$aic
}
Res$loglik <- loglik <- -opt$value
names(Res$pars) <- c("a","b","sd","rho")
Res$pars <- data.frame(t(Res$pars))
# Res$gamma <- gamma
ssb.tmp <- seq(from=0,to=max(ssb)*max.ssb.pred,length=100)
R.tmp <- sapply(1:length(ssb.tmp), function(i) SRF(ssb.tmp[i],a,b))
pred.data <- data.frame(SSB=ssb.tmp,R=R.tmp)
Res$pred <- pred.data
Res$k <- k <- length(opt$par)+1
Res$AIC <- -2*loglik+2*k
Res$AICc <- Res$AIC+2*k*(k+1)/(NN-k-1)
Res$BIC <- -2*loglik+k*log(NN)
return(Res)
}
### 西嶋加筆
# Allee effect (depensation)ありの再生産関係の推定用関数 (c.est=FALSEとすればfit.SRと同じ)
# 修正が必要
fit.SR2 <- function(SRdata,
SR="HS",
method="L2",
AR=1,
hessian=FALSE,
w=rep(1,length(SRdata$year)),
length=20, #parameter (a,b) の初期値を決めるときにgrid searchする数
c.est = TRUE #Allee effectを推定するかどうか(c>1でdepensation (Allee-like), c<1でcompensation)
){
argname <- ls()
arglist <- lapply(argname,function(xx) eval(parse(text=xx)))
names(arglist) <- argname
rec <- SRdata$R
ssb <- SRdata$SSB
N <- length(rec)
NN <- sum(w) #sample size for likelihood calculation
if (SR=="HS") SRF <- function(x,a,b,c) ifelse(x>b,b*a,a*b*(x/b)^c)
if (SR=="BH") SRF <- function(x,a,b,c) (a/b)/(1+1/(b*x)^c)
if (SR=="RI") SRF <- function(x,a,b,c) a/(b*exp(1))*(b*x)^c*exp(c*(1-b*x))
obj.f <- function(a,b,rho,c){
resid <- sapply(1:N,function(i) log(rec[i]) - log(SRF(ssb[i],a,b,c)))
resid2 <- NULL
for (i in 1:N) {
resid2[i] <- ifelse(i==1,resid[i], resid[i]-rho*resid[i-1])
}
if (method == "L2") {
sd <- sqrt(sum(w*resid2^2)/(NN-rho^2))
sd2 <- c(sd/sqrt(1-rho^2), rep(sd,N-1))
obj <- -sum(w*dnorm(resid2,0,sd2,log=TRUE))
} else {
sd <- sum(abs(w*resid2))/(NN-rho^2)
sd2 <- c(sd/sqrt(1-rho^2), rep(sd,N-1))
obj <- -sum(w*sapply(1:N, function(i){-log(2*sd2[i])-abs(resid2[i]/sd2[i])}))
}
return(obj)
}
a.range <- range(rec/ssb)
b.range <- range(1/ssb)
if (SR == "HS") b.range <- range(ssb)
grids <- as.matrix(expand.grid(
seq(a.range[1],a.range[2],len=length),
seq(b.range[1],b.range[2],len=length)
))
init <- as.numeric(grids[which.min(sapply(1:nrow(grids),function(i) obj.f(grids[i,1],grids[i,2],0,1))),])
init[1] <- log(init[1])
init[2] <- ifelse (SR == "HS",-log(max(0.000001,(max(ssb)-min(ssb))/max(init[2]-min(ssb),0.000001)-1)),log(init[2]))
if (AR != 0 || isTRUE(c.est)) init[3] <- 0
if (AR != 0 && isTRUE(c.est)) init[4] <- 0
if (SR == "HS") {
if (AR == 0) {
if (c.est) {
obj.f2 <- function(x) obj.f(exp(x[1]),min(ssb)+(max(ssb)-min(ssb))/(1+exp(-x[2])),0,exp(x[3]))
} else {
obj.f2 <- function(x) obj.f(exp(x[1]),min(ssb)+(max(ssb)-min(ssb))/(1+exp(-x[2])),0,1)
}
} else {
if (c.est) {
obj.f2 <- function(x) obj.f(exp(x[1]),min(ssb)+(max(ssb)-min(ssb))/(1+exp(-x[2])),1/(1+exp(-x[3])),exp(x[4]))
} else {
obj.f2 <- function(x) obj.f(exp(x[1]),min(ssb)+(max(ssb)-min(ssb))/(1+exp(-x[2])),1/(1+exp(-x[3])),1)
}
}
} else {
if (AR == 0) {
if (c.est) {
obj.f2 <- function(x) obj.f(exp(x[1]),exp(x[2]),0,exp(x[3]))
} else {
obj.f2 <- function(x) obj.f(exp(x[1]),exp(x[2]),0,1)
}
} else {
if (c.est) {
obj.f2 <- function(x) obj.f(exp(x[1]),exp(x[2]),1/(1+exp(-x[3])),exp(x[4]))
} else {
obj.f2 <- function(x) obj.f(exp(x[1]),exp(x[2]),1/(1+exp(-x[3])),1)
}
}
}
opt <- optim(init,obj.f2)
opt <- optim(opt$par,obj.f2,method="BFGS",hessian=hessian)
Res <- list()
Res$input <- arglist
Res$opt <- opt
a <- exp(opt$par[1])
b <- ifelse(SR=="HS",min(ssb)+(max(ssb)-min(ssb))/(1+exp(-opt$par[2])),exp(opt$par[2]))
rho <- ifelse(AR==0,0,1/(1+exp(-opt$par[3])))
c <- ifelse(c.est, exp(rev(opt$par)[1]),1)
resid <- sapply(1:N,function(i) log(rec[i]) - log(SRF(ssb[i],a,b,c)))
resid2 <- NULL
for (i in 1:N) {
resid2[i] <- ifelse(i == 1,resid[i], resid[i]-rho*resid[i-1])
}
sd <- ifelse(method=="L2",sqrt(sum(w*resid2^2)/(NN-rho^2)),sqrt(2)*sum(abs(w*resid2))/(NN-rho^2))
Res$resid <- resid
Res$resid2 <- resid2
Res$pars <- c(a,b,sd,rho,c)
if (method!="L2") {
if (AR!=0) {
arres <- ar(resid,aic=FALSE,order.max=1)
Res$pars[3] <- sqrt(arres$var.pred)
Res$pars[4] <- arres$ar
}
}
Res$loglik <- loglik <- -opt$value
names(Res$pars) <- c("a","b","sd","rho","c")
Res$pars <- data.frame(t(Res$pars))
ssb.tmp <- seq(from=0,to=max(ssb)*1.3,length=100)
R.tmp <- sapply(1:length(ssb.tmp), function(i) SRF(ssb.tmp[i],a,b,c))
pred.data <- data.frame(SSB=ssb.tmp,R=R.tmp)
Res$pred <- pred.data
Res$k <- k <- length(opt$par)+1
Res$AIC <- -2*loglik+2*k
Res$AICc <- Res$AIC+2*k*(k+1)/(NN-k-1)
Res$BIC <- -2*loglik+k*log(NN)
return(Res)
}
plot.kobe <- function(vpares,Bmsy,Umsy,Blim=NULL,Bban=NULL,plot.history=FALSE,is.plot=FALSE,pickU="",pickB="",
ylab.tmp=ifelse(yaxis=="U","U/Umsy","F/Fmsy"),
xlab.tmp="SSB/SSBmsy",title.tmp="",HCR=NULL,
yaxis="U" # y軸になにをとるか。Uの場合は漁獲率。Fの場合は F=-log(1-U)として計算したfishing mortality
){ # HCR=list(beta=0.8)
if (is.null(vpares$wcaa)) vpares$wcaa <- vpares$input$dat$caa * vpares$input$dat$waa
vpares$TC.MT <- as.numeric(colSums(vpares$wcaa))
U_history <- as.numeric(vpares$TC.MT)/as.numeric(colSums(vpares$baa,na.rm=T))
F_history <- -log(1-U_history)
F_msy <- -log(1-Umsy)
UBdata <- data.frame(years=as.numeric(colnames(vpares$baa)),
U=U_history/Umsy,
B=as.numeric(colSums(vpares$ssb))/Bmsy,
F=F_history/F_msy)
x <- UBdata$B
if(yaxis=="U") y <- UBdata$U else y <- UBdata$F
tmp <- x>0 & y>0
x <- x[tmp]
y <- y[tmp]
UBdata <- UBdata[tmp,]
if(!is.null(Blim)){
Blim.percent <- Blim/Bmsy
}
else{
Blim.percent <- 0.5
}
plot(x,
y,type="n",xlim=c(0,ifelse(max(x)<2,2,max(x,na.rm=T))),
ylim=c(0,ifelse(max(y,na.rm=T)<3,3,max(y,na.rm=T))),
cex=c(1,rep(1,length(y)-2),3),ylab=ylab.tmp,xlab=xlab.tmp)
polygon(c(-1,1,1,-1),c(-1,-1,1,1),col="khaki1",border=NA)
polygon(c(1,6,6,1),c(-1,-1,1,1),col="olivedrab2",border=NA)
polygon(c(1,6,6,1),c(1,1,6,6),col="khaki1",border=NA)
polygon(c(-1,Blim.percent,Blim.percent,-1),c(1,1,6,6),col="indianred1",border=NA)
polygon(c(Blim.percent,1,1,Blim.percent),c(1,1,6,6),col="tan1",border=NA)
polygon(c(-1,Blim.percent,Blim.percent,-1),c(-1,-1,1,1),col="khaki2",border=NA)
polygon(c(Blim.percent,1,1,Blim.percent),c(-1,-1,1,1),col="khaki1",border=NA)
axis(side=1:2)
if(!is.null(HCR)){
lines(c(Bban/Bmsy,Blim/Bmsy,6),c(0,0.8,0.8),lty=2)
}
# points(x,y,type="o",pch=c(3,rep(1,length(y)-2),20),col=c(1,rep(1,length(y)-2),1),cex=c(1,r
points(x,y,type="l",pch=20,col=1,lwd=1)
points(x,y,type="p",pch=20,col=gray(c(seq(from=0.7,to=0,length=length(x)))),cex=1.2)
points(rev(x)[1],rev(y)[1],type="p",pch=20,cex=2.5)
title(title.tmp,adj=0.8,line=-2)
if(isTRUE(plot.history)){
plot(UBdata$years,UBdata$U,type="b",ylab="U/Umsy",xlab="Year",ylim=c(0,max(y)))
abline(h=1)
plot(UBdata$years,UBdata$F,type="b",ylab="F/Fmsy",xlab="Year",ylim=c(0,max(y)))
abline(h=1)
plot(UBdata$years,UBdata$B,type="b",ylab="SSB/SSBmsy",xlab="Year",ylim=c(0,max(y)))
abline(h=1); abline(h=Blim.percent,lty=2)
}
invisible(UBdata)
}
plot.waa <- function(vres){
lm.list <- list()
nage <- nrow(vres$naa)
col.tmp <- rainbow(nage)
logx <- log(unlist(vres$naa))
logy <- log(unlist(vres$input$dat$waa))
ages <- as.numeric(rep(rownames(vres$naa),ncol(vres$naa)))
u.age <- unique(ages)
plot(logx,logy,col=col.tmp[1+ages],xlab="log(N)",ylab="log(weight)")
for(i in 1:length(u.age)){
tmp <- ages==u.age[i] & logy>-Inf & logx>-Inf
if(sum(tmp,na.rm=TRUE)>0){
lm.list[[i]] <- lm(logy[tmp]~logx[tmp])
l.type <- ifelse(summary(lm.list[[i]])$coeff[2,4]<0.05,1,2)
if(!is.na(l.type)) abline(lm.list[[i]],col=col.tmp[1+ages[i]],lty=l.type)
}
}
title(vres$stockid,line=0.2)
legend("bottomleft",lty=c(1:2,rep(1,nage)),
col=c(1,1,col.tmp),
legend=c("p<0.05","p>0.05",paste("Age",u.age)))
return(lm.list)
}
### parametric bootstrap usnig fit.SR
boot.SR <- function(Res,n=100,seed=1){
N <- length(Res$input$SRdata$year)
# if (Res$input$SR=="HS") SRF <- function(x,a,b,gamma=Res$gamma) a*(x+sqrt(b^2+gamma^2/4)-sqrt((x-b)^2+gamma^2/4))
if (Res$input$SR=="HS") SRF <- function(x,a,b) ifelse(x>b,b*a,x*a)
if (Res$input$SR=="BH") SRF <- function(x,a,b) a*x/(1+b*x)
if (Res$input$SR=="RI") SRF <- function(x,a,b) a*x*exp(-b*x)
sd <- sapply(1:N, function(i) ifelse(i==1,Res$pars$sd/sqrt(1-Res$pars$rho^2),Res$pars$sd))
set.seed(seed)
lapply(1:n, function(j){
N <- length(Res$input$SRdata$SSB)
resids <- rnorm(N,0,sd)
pred <- obs <- resid0 <- numeric(N)
ssb <- Res$input$SRdata$SSB
for(i in 1:N){
pred[i] <- SRF(ssb[i],Res$pars$a,Res$pars$b)
if (i==1) {
obs[i] <- pred[i]*exp(resids[i])
} else {
obs[i] <- pred[i]*exp(Res$pars$rho*resid0[i-1])*exp(resids[i])
}
resid0[i] <- log(obs[i]/pred[i])
}
res.b <- Res
res.b$input$SRdata$R <- obs
res.b <- do.call(fit.SR, res.b$input)
return(res.b)
})
}
### profile likelihood
prof.lik <- function(Res,a=Res$pars$a,b=Res$pars$b,sd=Res$pars$sd,rho=Res$pars$rho) {
SRdata <- Res$input$SRdata
rec <- SRdata$R
ssb <- SRdata$SSB
N <- length(rec)
SR <- Res$input$SR
gamma <- Res$gamma
method <- Res$input$method
w <- Res$input$w
# if (SR=="HS") SRF <- function(x,a,b) a*(x+sqrt(b^2+gamma^2/4)-sqrt((x-b)^2+gamma^2/4))
if (SR=="HS") SRF <- function(x,a,b) ifelse(x>b,b*a,x*a)
if (SR=="BH") SRF <- function(x,a,b) a*x/(1+b*x)
if (SR=="RI") SRF <- function(x,a,b) a*x*exp(-b*x)
resid <- sapply(1:N,function(i) log(rec[i]) - log(SRF(ssb[i],a,b)))
resid2 <- NULL
for (i in 1:N) {
resid2[i] <- ifelse(i==1,resid[i], resid[i]-rho*resid[i-1])
}
obj <- NULL
if (method == "L2") {
for (i in 1:N) {
if (i==1) {
obj <- c(obj,-0.5*log(2*pi)-log(sd^2/(1-rho^2))-resid2[i]^2/(2*sd^2/(1-rho^2)))
} else {
obj <- c(obj, -0.5*log(2*pi)-0.5*log(sd^2)-resid2[i]^2/(2*sd^2))
}
}
} else {
for (i in 1:N) {
if (i==1) {
obj <- c(obj,-log(2*sqrt(sd^2/(1-rho^2)))-abs(resid2[i])/sqrt(sd^2/(1-rho^2)))
} else {
obj <- c(obj, -log(2*sd)-abs(resid2[i])/sd)
}
}
}
obj <- sum(w*obj) # exact likelihood
return(exp(obj))
}
## ちょっと複雑なkobe.plot
# foutsが複数の将来予測の結果。brefsは複数の管理基準値
get.kobemat2 <- function(fouts,brefs,xlim=NULL,target.prob=0.5){
# brefs <- sort(brefs)
years <- as.numeric(rownames(fouts[[1]]$vssb))
probs <- matrix(0,length(years),length(fouts))
for(j in 1:ncol(brefs)){
probs <- probs + foreach(i=1:length(fouts),.combine=cbind) %do%
as.numeric(rowMeans(fouts[[i]]$vssb > brefs[i,j])>target.prob)
}
if(is.null(xlim)) xlim <- range(years)
plot(range(years),
range(0.5,nrow(brefs)+1.5),
type="n",xlab="Years",cex=3,xlim=xlim,
ylab="Strategies",yaxt="n")
abline(h=1:ncol(brefs),v=years,col="gray")
axis(side=2,at=1:nrow(brefs),label=rownames(brefs))
require(RColorBrewer)
cols <- brewer.pal(ncol(brefs), "Paired")
for(i in 1:length(fouts)){
points(years,rep(i,length(years)),
col=cols[probs[,i]],pch=20,cex=3)
}
legend("topright",pch=20,cex=1,col=cols,ncol=ceiling(ncol(brefs)/2),
legend=paste("Prob(B>",colnames(brefs),")>",round(target.prob*100),"%"))
}
Generation.Time <- function(vpares,
maa.year=2014:2015,
M.year=2014:2015,
Plus = 19
){
maa <- vpares$input$dat$maa
maa <- rowMeans(maa[,colnames(maa) %in% maa.year,drop=F],na.rm=T)
maa <- maa[!is.na(maa)]
M <- vpares$input$dat$M
M <- rowMeans(M[,colnames(M) %in% M.year,drop=F],na.rm=T)
M <- M[!is.na(M)]
age <- as.numeric(names(maa))
maa <- c(maa, rep(1,Plus))
M <- c(M, rep(M[length(M)],Plus))
age <- c(age, max(age)+1:Plus)
A <- length(M)
L <- c(1,exp(-cumsum(M[-A])))
G <- sum(age*L*maa)/sum(L*maa)
return(G)
}
###############################
#### 資源量の上積みグラフを書く
###############################
plotBfish <- function(tres0,vpares, # SR.estの結果
b.target,ssb.max=Inf,
biomass.scale=1000){
ssb <- tres0$ssb.mean/biomass.scale
tmp <- substr(colnames(tres0),1,5)=="TB-MA"
tb <- tres0[,tmp]/biomass.scale
tb2 <- sapply(1:ncol(tb),function(x) apply(tb[,1:x,drop=F],1,sum,na.rm=T))
tmp <- substr(colnames(tres0),1,5)=="TC-MA"
tc <- tres0[,tmp]/biomass.scale
tc2 <- sapply(1:ncol(tc),function(x) apply(tc[,1:x,drop=F],1,sum,na.rm=T))
# library(png)
# {if(file.exists(pngfile)) image <- readPNG(pngfile)
# else image <- NULL}
year.tmp <- rev(colnames(vpares$ssb))[1:5]
range1 <- c(0,min(ssb.max,max(ssb)))
range2 <- range(as.data.frame(vpares$ssb)[as.character(year.tmp)])
col.tmp1 <- rgb(40/255,96/255,163/255,seq(from=0.1,to=0.9,length=ncol(tc)))
col.tmp2 <- rgb(100/255,200/255,44/255,seq(from=0.1,to=0.9,length=ncol(tc)))
### plot of SSB
# tb3 <- tb2[which(ssb<ssb.max),]
matplot(ssb,tb2,type="n",ylab=paste("Total biomass (",biomass.scale," MT)",sep=""),xaxs="i",yaxs="i",
xlab="SSB (1000MT)", xlim=range1,ylim=c(0,max(tb2)))
# 過去の時系列
matpoints(ssb,tb2[,1],type="l",lwd=2,col="gray",lty=3)
# points(x <- colSums(vpares$ssb)/biomass.scale,
# y <- colSums(vpares$baa)/biomass.scale,type="o",
# col=gray(c(seq(from=0.7,to=0,length=length(x)))),pch=20,cex=1.2,
# lwd=3)
# text(x[1],y[1],colnames(x)[1],adj=0)
#text(rev(x[1]),rev(y)[1],rev(colnames(x))[1],adj=0)
text(rep(ssb[10],ncol(tb2)),
tb2[10,],
paste(0:(ncol(tb2)-1),"y/o"))#,round(apply(vpares$input$dat$waa,1,mean),0)," g"),cex=1)
ssb.hist <- range(colSums(vpares$ssb)/biomass.scale)
polygon(c(ssb.hist,rev(ssb.hist)),c(0,0,max(tb2)*10,max(tb2)*10),col=gray(0.8),border=NA)
abline(v=b.target,col="gray")
## 積み上げグラフ
non.na <- !is.na(ssb)
for(i in 1:ncol(tb2)) menplot(ssb[non.na], cbind(0,tb2)[non.na,i:(i+1)],col=col.tmp1[i],border=NA)
# title("Total biomass",line=-1,adj=0.1)
## catch
{if(!is.null(vpares$wcaa)) wcatch <- as.numeric(colSums(vpares$wcaa))
else{
wcatch <- as.numeric(colSums(vpares$input$dat$caa * vpares$input$dat$waa,na.rm=T))
}}
matplot(ssb,tc2,type="n",,xaxs="i",yaxs="i",ylab=paste("Catch (",biomass.scale," MT)",sep=""),
xlab="SSB (1000MT)",
ylim=c(0,max(tc2)*1.2),xlim=range1)
polygon(c(ssb.hist,rev(ssb.hist)),c(0,0,max(tc2)*10,max(tc2)*10),col=gray(0.8),border=NA)
abline(v=b.target,col="gray")
for(i in 1:ncol(tc2)) menplot(ssb[non.na], cbind(0,tc2)[non.na,i:(i+1)],col=col.tmp2[i],border=NA)
text(rep(ssb[10],ncol(tc2)),
tc2[10,],
paste(0:(ncol(tc2)-1),"y/o"),col="darkgreen")#,round(apply(vpares$input$dat$waa,1,mean),0)," g"),cex=1)
# title("Total catch",line=-1,adj=0.1)
if(0){
plot(ssb,tres0$fmulti,type="n",lwd=2,xlim=range1,ylab="Fishing efforts",xlab="SSB (1000MT)")
polygon(c(ssb.hist,rev(ssb.hist)),c(0,0,10,10),col=gray(0.8),border=NA)
abline(v=b.target,col="gray")
points(ssb,tres0$fmulti,type="l",lwd=2,xlim=range1,ylab="Fishing efforts")
par(new=T)
matplot(ssb,cbind(tres0$ssb.CV,tres0$catch.CV),xlim=range1,ylim=c(0,2),type="l",lty="22",col=3:4,lwd=2,axes=F,ylab="",xlab="")
axis(side=4)
mtext(side=4,"CV",cex=1,line=2.3)
}
}
############## 岡村さん作成関数 ##################3
future.vpa1 <- function(
vpares,
multi=1,
nyear=50, # 将来予測の年数
ABC.year=2018, # ABCを計算する年
waa.year=2015:2016, # 生物パラメータの参照年
maa.year=2015:2016,
M.year=2015:2016,
seed=1,
N=100,
naa0=NULL,
eaa0=NULL,
beta=1,
delta=0,
Blim=0,
Bban=0,
Pope=FALSE,
ssb0=NULL,
faa0=NULL,
# recfuncに対する引数
rec.arg=list(a=res1$pars[1],b=res1$pars[2],gamma=res1$gamma,
sd=res1$pars[3],bias.correction=TRUE,rho=res1$pars[4],resid=res1$resid)
){
# print(multi)
argname <- ls()
arglist <- lapply(argname,function(x) eval(parse(text=x)))
names(arglist) <- argname
set.seed(seed)
lag <- as.numeric(rownames(vpares$input$dat$caa))[1]
nY <- ncol(vpares$input$dat$caa)
final.year <- as.numeric(colnames(vpares$input$dat$caa)[nY])
waa <- vpares$input$dat$waa
waa <- rowMeans(waa[,colnames(waa) %in% waa.year])
maa <- vpares$input$dat$maa
maa <- rowMeans(maa[,colnames(maa) %in% maa.year])
M <- vpares$input$dat$M
A <- nrow(M)
M <- rowMeans(M[,colnames(M) %in% M.year])
waa <- array(waa,dim=c(A,nyear+1,N))
maa <- array(maa,dim=c(A,nyear+1,N))
M <- array(M,dim=c(A,nyear+1,N))
naa <- baa <- ssb <- faa <- caa <- wcaa <- array(NA,dim=c(A,nyear+1,N))
if (is.null(faa0)) Fc.at.age <- vpares$Fc.at.age else Fc.at.age <- faa0
# 1st year
if (is.null(naa0)) {
waa[,1,] <- vpares$input$dat$waa[,nY]
maa[,1,] <- vpares$input$dat$maa[,nY]
M[,1,] <- vpares$input$dat$M[,nY]
naa[,1,] <- vpares$naa[,nY]
faa[,1,] <- vpares$faa[,nY]
baa[,1,] <- naa[,1,]*waa[,1,]
ssb[,1,] <- baa[,1,]*maa[,1,]
faa[,2,] <- Fc.at.age
faa[,-(1:2),] <- multi*Fc.at.age
} else {
naa[,1,] <- naa0
faa[,,] <- multi*Fc.at.age
baa[,1,] <- naa[,1,]*waa[,1,]
ssb[,1,] <- baa[,1,]*maa[,1,]
}
eaa <- matrix(NA,nrow=nyear+1,ncol=N)
sd2 <- sqrt(rec.arg$sd^2/(1-rec.arg$rho^2))
if (is.null(eaa0)) eaa[1,] <- rec.arg$resid[length(rec.arg$resid)] else eaa[1,] <- eaa0
eaa[-1,] <- rnorm(nyear*N,0,rec.arg$sd)
if(class(ssb0)=="matrix") ssb0 <- array(ssb0,dim=c(A,1,N))
if (is.null(ssb0) & (2 - lag <= 0)) {
ssb0 <- array(NA, dim=c(A,abs(2-lag)+1,N))
for (j in 1:(abs(2-lag)+1)){
ssb0[,j,] <- vpares$ssb[,nY-j]
}
}
alpha <- array(1,dim=c(A,nyear+1,N))
# 2nd year and onward
for (i in 2:(nyear+1)){
eaa[i,] <- rec.arg$rho*eaa[i-1,]+eaa[i,]
naa[2:A,i,] <- naa[1:(A-1),i-1,]*exp(-M[1:(A-1),i-1,]-alpha[1:(A-1),i-1,]*faa[1:(A-1),i-1,])
naa[A,i,] <- naa[A,i,]+naa[A,i-1,]*exp(-M[A,i-1,]-alpha[A,i-1,]*faa[A,i-1,])
ssb[,i,] <- naa[,i,]*waa[,i,]*maa[,i,]
if(i-lag > 0) SSB <- colSums(ssb[,i-lag,,drop=FALSE],na.rm=TRUE) else SSB <- colSums(ssb0[,i-1,,drop=FALSE])
naa[1,i,] <- HS(SSB,rec.arg$a,rec.arg$b,rec.arg$gamma,HStype="HS")*exp(eaa[i,]-0.5*sd2^2)
baa[,i,] <- naa[,i,]*waa[,i,]
ssb[,i,] <- baa[,i,]*maa[,i,]
Bcur <- colSums(ssb[,i,,drop=FALSE],na.rm=TRUE)
alpha[,i,] <- beta*matrix(ifelse(Bcur > Blim, 1, ifelse(Bcur > Bban, ((Bcur-Bban)/(Blim-Bban))^delta, 0)),byrow=TRUE,nrow=A,ncol=N)
}
if (Pope) caa <- (1-exp(-alpha*faa))*exp(-M/2)*naa else caa <- naa*(1-exp(-alpha*faa-M))*alpha*faa/(alpha*faa+M)
wcaa <- caa*waa
vwcaa <- apply(wcaa,c(2,3),sum,na.rm=T)
res <- list(beta=beta,delta=delta,alpha=alpha,Blim=Blim,Bban=Bban,waa=waa,maa=maa,M=M,naa=naa,baa=baa,ssb=ssb,faa=faa,caa=caa,wcaa=wcaa,vwcaa=vwcaa,eaa=eaa,multi=multi,input=arglist)
}
#
HS <- function(x,a,b,gamma1=0.001,HStype="HS") if (HStype=="Mesnil") a*(x+sqrt(b^2+(gamma1^2)/4)-sqrt((x-b)^2+(gamma1^2)/4)) else ifelse(x > b, a*b, a*x)
##
est.MSY2 <- function(vpares,N=1000,res1=NULL,sim0=NULL,nyear=NULL,pgy=0.9,lim=0.6,ban=0.1,mY=5,long.term=20,
Fmsy.max=3, # current FがFmsyに比べて小さすぎる場合、うまく収束しない場合があります。そのときはこのオプションでFmsy.max=10とかしてください。
Fmsy.step=0.1,thin=1,inc=1,SRtype="L2",fm=5,tol=NULL,
AutoCor=FALSE,# 関数内部で自己相関係数を推定するか "future.vpa"を使う場合はどちらでも良い
AutoCorOut=FALSE, # フィットさせたあと残差の自己相関を計算する場合
current.resid=0, # 最近年何年分の自己相関を平均するか
future.function.name="future.vpa1",seed=1){
if (is.null(tol)) tol <- .Machine$double.eps^0.25
Ccur <- sum(tail(t(vpares$input$dat$caa*vpares$input$dat$waa),1),na.rm=TRUE)
Blim.cur <- vpares$Blim
A <- nrow(vpares$input$dat$caa)
if (is.na(vpares$Fc.at.age[length(vpares$Fc.at.age)])){
vpares$input$dat$caa <- vpares$input$dat$caa[-A,]
vpares$input$dat$waa <- vpares$input$dat$waa[-A,]
vpares$input$dat$maa <- vpares$input$dat$maa[-A,]
vpares$input$dat$M <- vpares$input$dat$M[-A,]
vpares$naa <- vpares$naa[-A,]
vpares$faa <- vpares$faa[-A,]
vpares$Fc.at.age <- vpares$Fc.at.age[-A]
vpares$ssb[is.na(vpares$ssb)] <- 0
A <- A-1
}
vpares$Fc.at.age <- fm*vpares$Fc.at.age
SRdata <- get.SRdata(vpares)
# fit SR
if (is.null(res1) && is.null(sim0)){
res0 <- estSR(SRdata,SR="HS",type=SRtype,Length=20,rho.range=0,AutoCor=FALSE)
if (AutoCor){
res1 <- estSR(SRdata,SR="HS",type=SRtype,Length=20,rho.range=0,AutoCor=AutoCor)
if (res0$AICc <= res1$AICc) res1 <- res0
# if(!is.null(res0$aic)) if(res0$aic <= res1$aic) res1 <- res0
# if(!is.null(res0$AICc)) if(res0$AICc <= res1$AICc) res1 <- res0
} else{
res1 <- res0
}
}
if (class(res1$pars)!="numeric") res1$pars <- as.numeric(res1$pars)
if (AutoCorOut){
ar1 <- ar(ts(res1$resid),order.max=1)
rho <- ar1$ar
if (length(rho)==0) rho <- 0
if (abs(rho) >= 0.99) rho <- sign(rho)*0.99
res1$pars[3] <- sqrt(ar1$var.pred)
res1$pars[4] <- rho
}
if (current.resid > 0) w.recent <- mean(rev(res1$resid)[1:current.resid]) else w.recent <- 0
# Initial Setting
lag <- as.numeric(rownames(vpares$input$dat$caa))[1]
Pope <- vpares$input$Pope
future.vpa1 <- get(future.function.name)
years <- sort(as.numeric(rev(names(vpares$naa))[1:5]))
Surv <- exp(-tail(t(vpares$input$dat$M),1))
L <- cumprod(Surv)
L[A-1] <- L[A-1]/(1-Surv[A])
L <- c(1,L[1:(A-1)])
GT <- Generation.Time(vpares,maa.year=years,M.year=years) # Generation Time
GT2 <- round(GT*2)
det.naa0 <- res1$pars[1]*res1$pars[2]*L
waa <- vpares$input$dat$waa
waa <- rowMeans(waa[,colnames(waa) %in% years])
maa <- vpares$input$dat$maa
maa <- rowMeans(maa[,colnames(maa) %in% years])
det.B0 <- sum(det.naa0*waa*maa)
if(is.null(nyear)) nyear <- round(GT*long.term)
if (is.null(sim0)){
sim0 <- future.vpa1(vpares,
multi=0,
nyear=nyear, # 将来予測の年数
N=N, # 確率的計算の繰り返し回数
ABC.year=max(years)+1, # ABCを計算する年
waa.year=years, # 生物パラメータの参照年
maa.year=years,
M.year=years,
seed=seed,
naa0=det.naa0,
Pope=Pope,
# recfuncに対する引数
rec.arg=list(a=res1$pars[1],b=res1$pars[2],gamma=res1$gamma,sd=res1$pars[3],bias.correction=TRUE,rho=res1$pars[4],resid=res1$resid)
)
sim1 <- future.vpa1(vpares,
multi=1,
nyear=nyear, # 将来予測の年数
N=1, # 確率的計算の繰り返し回数
ABC.year=max(years)+1, # ABCを計算する年
waa.year=years, # 生物パラメータの参照年
maa.year=years,
M.year=years,
seed=seed,
naa0=det.naa0,
Pope=Pope,
# recfuncに対する引数
rec.arg=list(a=res1$pars[1],b=res1$pars[2],gamma=res1$gamma,sd=res1$pars[3],bias.correction=TRUE,rho=res1$pars[4],resid=res1$resid)
)
} else{
farg <- sim0$input
farg$N <- N
farg$nyear <- nyear
farg$multi <- 0
farg$ABC.year <- max(years)+1
farg$naa0 <- det.naa0
if(!is.null(farg$pre.catch)){
farg$pre.catch <- NULL # pre.catchオプションがあるとうまくいかないのでなかったことにする
cat("notice: option \"pre.catch\" is turned off in estimating MSY.\n")
}
if(!is.null(farg$rec.new)){
farg$rec.new <- NULL # rec.newプションがあるとうまくいかないのでなかったことにする
cat("notice: option \"rec.new\" is turned off in estimating MSY.\n")
}
farg$add.year <- 1
farg$is.plot <- FALSE
farg$silent <- TRUE
farg$det.run <- FALSE
sim0 <- do.call(future.vpa1,farg)
farg$N <- 2
farg$multi <- 1
sim1 <- do.call(future.vpa1,farg)
}
## MSY推定
farg <- sim1$input
nY <- nyear+1
eyear <- mY+(lag > 0)*(lag-1)
syfunc <- function(x,farg,nyear=50,N=100,eyear=4,naa0=NULL,eaa0=NULL,ssb0=NULL,faa0=NULL,sd=NULL){
farg$multi <- x
farg$N <- N
farg$nyear <- nyear
farg$naa0 <- naa0
farg$eaa0 <- eaa0
farg$ssb0 <- ssb0
farg$faa0 <- faa0
if (!is.null(sd)) farg$rec.arg$sd <- sd
fout <- do.call(future.vpa1,farg)
nY <- nyear+1
out <- list(catch=fout$vwcaa[(nY-(eyear-1)):nY,,drop=FALSE],ssb=fout$ssb[,(nY-(eyear-1)):nY,,drop=FALSE],naa=fout$naa[,(nY-(eyear-1)):nY,,drop=FALSE],baa=fout$baa[,(nY-(eyear-1)):nY,,drop=FALSE],eaa=fout$eaa[(nY-(eyear-1)):nY,,drop=FALSE])
return(out)
}
F.multi <- seq(0,Fmsy.max,by=Fmsy.step)
N0 <- sim0$naa[,nY,]
e0 <- sim0$eaa[nY,]
if(lag==0) SSB0 <- NULL else SSB0 <- sim0$ssb[,nY-(lag-1),]
FSYest <- lapply(F.multi, function(x) syfunc(x,farg,N=round(N/thin),nyear=nyear,eyear=eyear,naa0=N0[,1:round(N/thin)],eaa=e0[1:round(N/thin)],ssb0=SSB0[,1:round(N/thin)]))
FSYest.c <- sapply(1:length(F.multi), function(i) mean(FSYest[[i]]$catch))
num.msy <- which.max(FSYest.c)
num.msy0 <- num.msy
Fmsy.multi <- F.multi[num.msy]
obj.msy <- function(x) -mean(syfunc(x,farg,N=N,nyear=nyear,eyear=eyear,naa0=N0,eaa0=e0,ssb0=SSB0)$catch)
res.msy <- optimize(obj.msy, pmin(pmax(c(Fmsy.multi-inc*Fmsy.step,Fmsy.multi+inc*Fmsy.step),min(F.multi)),max(F.multi)),tol=tol)
Fmsy.multi <- res.msy$minimum
MSY <- -res.msy$objective
MSYres <- syfunc(Fmsy.multi,farg,N=N,nyear=nyear,eyear=eyear,naa0=N0,eaa0=e0,ssb0=SSB0)
ssb.msy <- mean(apply(MSYres$ssb,c(2,3),sum,na.rm=TRUE))
## PGY推定
id.pgy0 <- num.msy0:length(F.multi)
id.pgy <- which.min((FSYest.c[id.pgy0] - pgy*MSY)^2)
pgy.low <- id.pgy0[id.pgy]
Flow <- F.multi[pgy.low]
N.m <- MSYres$naa[,eyear,]
e.m <- MSYres$eaa[eyear,]
if(lag==0) SSB.m <- NULL else SSB.m <- MSYres$ssb[,eyear-(lag-1),]
obj.pgy <- function(x) (mean(syfunc(x,farg,N=N,nyear=nyear,eyear=eyear,naa0=N.m,eaa0=e.m,ssb0=SSB.m)$catch)-pgy*MSY)^2
res.pgy <- optimize(obj.pgy, pmin(pmax(c(Flow-inc*Fmsy.step,Flow+inc*Fmsy.step),Fmsy.multi),max(F.multi)),tol=tol)
Flow.multi <- res.pgy$minimum
PGYlow.res <- syfunc(Flow.multi,farg,N=N,nyear=nyear,eyear=eyear,naa0=N.m,eaa0=e.m,ssb0=SSB.m)
PGYlow <- mean(PGYlow.res$catch)
ssb.low <- mean(apply(PGYlow.res$ssb,c(2,3),sum,na.rm=T))
id.pgy0 <- 1:num.msy0
id.pgy <- which.min((FSYest.c[id.pgy0] - pgy*MSY)^2)
pgy.high <- id.pgy0[id.pgy]
Fhigh <- F.multi[pgy.high]
res.pgy <- optimize(obj.pgy, pmin(pmax(c(Fhigh-inc*Fmsy.step,Fhigh+inc*Fmsy.step),min(F.multi)),Fmsy.multi),tol=tol)
Fhigh.multi <- res.pgy$minimum
PGYhigh.res <- syfunc(Fhigh.multi,farg,N=N,nyear=nyear,eyear=eyear,naa0=N.m,eaa0=e.m,ssb0=SSB.m)
PGYhigh <- mean(PGYhigh.res$catch)
ssb.high <- mean(apply(PGYhigh.res$ssb,c(2,3),sum,na.rm=T))
## Bhs推定
if(res1$input$SR=="HS"){
det.Bhs <- res1$pars[2]
if (!is.null(SSB0)) SSB0.HS <- as.matrix(rowMeans(SSB0)) else SSB0.HS <- NULL
obj.HS <- function(x) sum(syfunc(x,farg,N=2,nyear=nyear,eyear=eyear,naa0=as.matrix(rowMeans(N0)),eaa0=NULL,ssb0=SSB0.HS,sd=0)$ssb[,eyear,1])-det.Bhs
res.HS <- uniroot(obj.HS, c(0,2*max(F.multi)),tol=tol)
FHS.multi <- res.HS$root
HSres <- syfunc(FHS.multi,farg,N=N,nyear=nyear,eyear=eyear,naa0=N0,eaa0=e0,ssb0=SSB0)
HScat <- mean(HSres$catch)
ssb.hs <- mean(apply(HSres$ssb,c(2,3),sum,na.rm=T))
}
else{
det.Bhs <- SSB0.HS <- obj.HS <- res.HS <- FHS.multi <- HSres <- HScat <- ssb.hs <- NULL
}
## target function
# 実際には最後の再生産関係の残差も入れてやる必要がある(自己相関を入れるときに必要)
## target.func <- function(x,farg,naa0=NULL,eaa0=NULL,ssb0=NULL,faa0=NULL,mY=5,N=1,seed=1,eyear=4,p=1,beta=1,delta=0,Blim=0,Bban=0,sd0=NULL){
target.func <- function(x,farg,naa0=NULL,eaa0=NULL,ssb0=NULL,faa0=NULL,mY=5,N=2,seed=1,eyear=4,p=1,beta=1,delta=0,Blim=0,Bban=0,sd0=NULL){
farg$multi <- x
farg$seed <- seed
farg$N <- N
farg$nyear <- mY
farg$naa0 <- p*naa0
farg$eaa0 <- eaa0
farg$ssb0 <- p*ssb0
farg$faa0 <- faa0
farg$beta <- beta
farg$delta <- delta
farg$Blim <- Blim
farg$Bban <- Bban
if(!is.null(farg$ABC.year)) farg$ABC.year <- farg$start.year
if (!is.null(sd0)) farg$rec.arg$sd <- sd0
fout <- do.call(future.vpa1,farg)
nY <- mY+1
out <- list(catch=fout$vwcaa[(nY-(eyear-1)):nY,,drop=FALSE],ssb=fout$ssb[,(nY-(eyear-1)):nY,,drop=FALSE],naa=fout$naa[,(nY-(eyear-1)):nY,,drop=FALSE],baa=fout$baa[,(nY-(eyear-1)):nY,,drop=FALSE],eaa=fout$eaa[(nY-(eyear-1)):nY,,drop=FALSE])
return(out)
}
## Blim0推定
Lim0.res <- LIMtoLOW <- list()
Flim.multi <- Lim0 <- ssb.lim0 <- PRT.lim <- numeric(length(lim))
PRT.range <- 1:(4*GT2)
N.m <- MSYres$naa[,eyear,]
e.m <- MSYres$eaa[eyear,]
if(lag==0) SSB.m <- NULL else SSB.m <- MSYres$ssb[,eyear-(lag-1),]
for (j in 1:length(lim)){
id.lim0 <- num.msy0:length(F.multi)
id.lim <- which.min((FSYest.c[id.lim0] - lim[j]*MSY)^2)
lim.num <- id.lim0[id.lim]
Flim <- F.multi[lim.num]
obj.lim <- function(x) (mean(syfunc(x,farg,N=N,nyear=nyear,eyear=eyear,naa0=N.m,eaa0=e.m,ssb0=SSB.m)$catch)-lim[j]*MSY)^2
res.lim <- optimize(obj.lim, pmin(pmax(c(Flim-inc*Fmsy.step,Flim+inc*Fmsy.step),Fmsy.multi),max(F.multi)),tol=tol)
Flim.multi[j] <- res.lim$minimum
Lim0.res[[j]] <- syfunc(Flim.multi[[j]],farg,N=N,nyear=nyear,eyear=eyear,naa0=N.m,eaa0=e.m,ssb0=SSB.m)
Lim0[j] <- mean(Lim0.res[[j]]$catch)
ssb.lim0[j] <- mean(apply(Lim0.res[[j]]$ssb,c(2,3),sum,na.rm=T))
## PRT
N.p <- Lim0.res[[j]]$naa[,eyear,]
e.p <- Lim0.res[[j]]$eaa[eyear,]
if(lag==0) SSB.p <- NULL else SSB.p <- Lim0.res[[j]]$ssb[,eyear-(lag-1),]
LIMtoLOW[[j]] <- target.func(Fmsy.multi,farg,mY=4*GT2,seed=seed,N=N,eyear=4*GT2,naa0=N.p,eaa0=e.p,ssb0=SSB.p)
LIMtoLOW.ssb <- sapply(PRT.range, function(x) mean(colSums(LIMtoLOW[[j]]$ssb[,x,])))
PRT.lim[j] <- min(which(LIMtoLOW.ssb >= ssb.low))
}
## PRT.lim <= mYを満たすものが一個もなくwarningを返すことがあるが、放置しています(使っていない計算結果なので)
nlim.est <- min(which(PRT.lim <= mY))
if(is.na(nlim.est) | nlim.est == Inf) nlim.est <- length(lim)
lim1 <- lim[nlim.est]
Lim1.res <- Lim0.res[[nlim.est]]
Flim1.multi <- Flim.multi[nlim.est]
Lim1 <- Lim0[nlim.est]
ssb.lim1 <- ssb.lim0[nlim.est]
PRT.lim1 <- PRT.lim[nlim.est]
## Bban0推定
Ban0.res <- BANtoLIM <- list()
Fban.multi <- Ban0 <- ssb.ban0 <- PRT.ban <- numeric(length(ban))
for (j in 1:length(ban)){
id.ban0 <- num.msy0:length(F.multi)
id.ban <- which.min((FSYest.c[id.ban0] - ban[j]*MSY)^2)
ban.num <- id.ban0[id.ban]
Fban <- F.multi[ban.num]
obj.ban <- function(x) (mean(syfunc(x,farg,N=N,nyear=nyear,eyear=eyear,naa0=N.m,eaa0=e.m,ssb0=SSB.m)$catch)-ban[j]*MSY)^2
res.ban <- optimize(obj.ban, pmin(pmax(c(Fban-inc*Fmsy.step,Fban+inc*Fmsy.step),Fmsy.multi),max(F.multi)),tol=tol)
Fban.multi[j] <- res.ban$minimum
Ban0.res[[j]] <- syfunc(Fban.multi[j],farg,N=N,nyear=nyear,eyear=eyear,naa0=N.m,eaa0=e.m,ssb0=SSB.m)
Ban0[j] <- mean(Ban0.res[[j]]$catch)
ssb.ban0[j] <- mean(apply(Ban0.res[[j]]$ssb,c(2,3),sum,na.rm=T))
## PRT
N.p <- Ban0.res[[j]]$naa[,eyear,]
e.p <- Ban0.res[[j]]$eaa[eyear,]
if(lag==0) SSB.p <- NULL else SSB.p <- Ban0.res[[j]]$ssb[,eyear-(lag-1),]
BANtoLIM[[j]] <- target.func(Fmsy.multi,farg,mY=4*GT2,seed=seed,N=N,eyear=4*GT2,naa0=N.p,eaa0=e.p,ssb0=SSB.p,delta=1,Blim=ssb.lim1,Bban=ssb.ban0[j])
BANtoLIM.ssb <- sapply(PRT.range, function(x) mean(colSums(BANtoLIM[[j]]$ssb[,x,])))
PRT.ban[j] <- min(which(BANtoLIM.ssb >= ssb.lim1))
}
nban.est <- min(which(PRT.ban <= mY))
if(is.na(nban.est) | nban.est == Inf) nban.est <- length(ban)
ban1 <- ban[nban.est]
Ban1.res <- Ban0.res[[nban.est]]
Fban1.multi <- Fban.multi[nban.est]
Ban1 <- Ban0[nban.est]
ssb.ban1 <- ssb.ban0[nban.est]
PRT.ban1 <- PRT.ban[nban.est]
## Btarget推定
## Btarget
Ftar.multi <- Fmsy.multi
eyear <- mY+(lag > 0)*(lag-1)
N.m <- MSYres$naa[,eyear,]
e.m <- MSYres$eaa[eyear,]
if(lag==0) SSB.m <- NULL else SSB.m <- MSYres$ssb[,eyear-(lag-1),]
TARres <- target.func(Ftar.multi,farg,mY=mY,seed=seed,N=N,eyear=mY,naa0=N.m,eaa0=e.m+w.recent,ssb0=SSB.m)
Btar <- mean(colSums(TARres$ssb[,mY,]))
# Blow 推定
N.low <- PGYlow.res$naa[,1+(lag>0)*(lag-1),]
e.low <- PGYlow.res$eaa[1+(lag>0)*(lag-1),]
# if(lag==0) SSB.low <- NULL else SSB.low <- PGYlow.res$ssb[,1,]
if(lag==0) SSB.low <- NULL else SSB.low <- PGYlow.res$ssb[,eyear-(lag-1),]
LOWres <- target.func(Flow.multi,farg,mY=mY,seed=seed,N=N,eyear=mY,naa0=N.low,eaa0=e.low+w.recent,ssb0=SSB.low)
Blow <- mean(colSums(LOWres$ssb[,mY,]))
P.low <- Blow/Btar
# Blim 推定
N.lim <- Lim1.res$naa[,1+(lag>0)*(lag-1),]
e.lim <- Lim1.res$eaa[1+(lag>0)*(lag-1),]
# if(lag==0) SSB.lim <- NULL else SSB.lim <- Lim1.res$ssb[,1,]
if(lag==0) SSB.lim <- NULL else SSB.lim <- Lim1.res$ssb[,eyear-(lag-1),]
LIMres <- target.func(Flim1.multi,farg,mY=mY,seed=seed,N=N,eyear=mY,naa0=N.lim,eaa0=e.lim+w.recent,ssb0=SSB.lim)
Blim <- mean(colSums(LIMres$ssb[,mY,]))
P.lim <- Blim/Btar
# Bban 推定
N.ban <- Ban1.res$naa[,1+(lag>0)*(lag-1),]
e.ban <- Ban1.res$eaa[1+(lag>0)*(lag-1),]
## if(lag==0) SSB.ban <- NULL else SSB.ban <- Ban1.res$ssb[,1,]
if(lag==0) SSB.ban <- NULL else SSB.ban <- Ban1.res$ssb[,eyear-(lag-1),]
BANres <- target.func(Fban1.multi,farg,mY=mY,seed=seed,N=N,eyear=mY,naa0=N.ban,eaa0=e.ban+w.recent,ssb0=SSB.ban)
Bban <- mean(colSums(BANres$ssb[,mY,]))
P.ban <- Bban/Btar
Pref <- c(P.low, P.lim, P.ban)
names(Pref) <- c("Low","Lim","Ban")
if(0){
# plotはしない
x.range <- range(SRdata$SSB,Btar,Bban)
plot(SRdata$SSB, SRdata$R,xlab="SSB",ylab="R",pch=16,col="blue",cex=1.2,xlim=x.range)
x.SSB <- seq(0,x.range[2],len=100)
lines(x.SSB,HS(x.SSB,res1$pars[1],res1$pars[2]),col="pink",lwd=3)
abline(v=Btar,col="green",lwd=3,lty=2)
abline(v=Blim,col="orange",lwd=3,lty=2)
abline(v=Bban,col="red",lwd=3,lty=2)
}
if(0){
argname <- ls()
arglist <- lapply(argname,function(x) eval(parse(text=x)))
names(arglist) <- argname
arglist$MSYres <- NULL
arglist$PGYlow.res <- NULL
invisible(arglist)
}
out <- list(
stockid=vpares$stockid,
seed=seed,
SRdata=SRdata,
res1=res1,
det.B0=det.B0,
det.Bhs=det.Bhs,
B0=mean(colSums(sim0$ssb[,nY,])),
Bmsy=ssb.msy,
Bpgy.low=ssb.low,
Bpgy.high=ssb.high,
Blim1=ssb.lim1,
Bban1=ssb.ban1,
Btar=Btar,
Blow=Blow,
Blim=Blim,
Bban=Bban,
Bhs=ssb.hs,
pgy=pgy,
lim=lim,
ban=ban,
fm=fm,
lag=lag,
eyear=eyear,
Fmsy=Fmsy.multi,
Flow=Flow.multi,
Fhigh=Fhigh.multi,
Flim=Flim1.multi,
Fban=Fban1.multi,
Fhs=FHS.multi,
farg=farg,
N=N,
GT=GT,
sim0=sim0,
w.recent=w.recent,
nyear=nyear,
MSYres=MSYres,
lim.est=lim1,
ban.est=ban1,
PRT.lim0=PRT.lim,
PRT.ban0=PRT.ban,
PRT.lim=PRT.lim1,
PRT.ban=PRT.ban1,
PGYlow.res=PGYlow.res,
PGYhigh.res=PGYhigh.res,
TARres=TARres,
LOWres=LOWres,
Lim0.res=Lim0.res,
Ban0.res=Ban0.res,
Lim1.res=Lim1.res,
Ban1.res=Ban1.res,
LIMres=LIMres,
BANres=BANres,
HSres=HSres,
Pref=Pref,
syfunc=syfunc,
target.func=target.func,
future.function.name=future.function.name,
Blim.cur=Blim.cur,
Ccur=Ccur
)
b.table <- unlist(out[c("Bmsy","Btar","Fmsy","Bpgy.low","Blow","Flow",
"Blim1","Blim","Flim","Bban1","Bban","Fban")])
# names(b.table) <- c("Bmsy\n(Equiribrium)","Bmsy\n(with AR)=Btarget","Fmsy",
# "Bpgy90%-low\n(Equiribrium)","Bpgy90%-low\n(with AR)=Blow","Flow")
# b.limit <- c(out$Blim0,out$Blim,out$Flim,out$Bban0,out$Bban,out$Fban)
b.table <- t(matrix(b.table,3,4))
# b.table[,1:2] <- b.table[,1:2]/1000
b.table <- cbind(apply(b.table[,1:2],2, function(x) round(as.numeric(x),0)),
round(b.table[,3],2))
w.recent2 <- ifelse(is.null(sim0$input$rec.arg$rho),NA,w.recent)
b.table <- rbind(b.table,c(NA,w.recent2,NA))
colnames(b.table) <- c("Equiribrium","with AR","Fref/Fcurrent")
rownames(b.table) <- c("Bmsy","B_pgy_90%_L","B_limit (B_pgy_60%_L)","B_ban (B_pgy_10%_L)","Recent residual")
out$summary <- b.table
return(out)
}
##
## ABC Calculation
##
### usage
## calc.beta(MSY.HS$input$msy,Ftar=refs$Fmsy,Btar=refs$Bmsy,Blim=refs$Blim,Bban=refs$Bban,N=1000)
calc.beta <- function(msy.input,Ftar=NULL,Btar=NULL,Blim=NULL,Bban=NULL,Blim.prob=0.9,Btar.prob=0.5,N=1000){
input.beta <- msy.input
input.beta$N <- N
input.beta$multi <- Ftar
## Blim, BbanをもとにしたHCRを使って将来予測を実施するためのオプションを追加
input.beta$HCR <- list(Blim=Blim,
Bban=Bban,
beta=1) # そのときのベータを探索するが、ここではとりあえず1としておく
input.beta$is.plot <- FALSE
input.beta$Frec <- list(stochastic=TRUE,
future.year=NULL, # NULLにしておくと将来予測の最終年の確率を見る
Blimit=Blim,
scenario="blimit", # 将来の親魚資源量をBlimitで指定した値を参照して決める
target.probs=100-Blim.prob*100) # Blimit「以下」になる確率を設定
fres.beta1 <- do.call(future.vpa,input.beta)
input.beta$Frec <- list(stochastic=TRUE,
future.year=NULL, # NULLにしておくと将来予測の最終年と判断する
Blimit=Btar,
scenario="blimit",target.probs=100-Btar.prob*100)
fres.beta2 <- do.call(future.vpa,input.beta)
beta <- min(fres.beta1$multi/Ftar, fres.beta2$multi/Ftar)
input.beta$multi <- beta*Ftar
input.beta$Frec <- NULL
fout <- do.call(future.vpa,input.beta)
cat("beta:",beta,"\n",
"Year:",rev(dimnames(fout$vssb)[[1]])[1],"\n",
"Prob(SSB>Btar):",mean(fout$vssb[nrow(fout$vssb),]>Btar)*100,"\n",
"Prob(SSB>Blim):",mean(fout$vssb[nrow(fout$vssb),]>Blim)*100,"\n")
return(beta)
}
## 岡村さん作成バージョン。新しい関数に差し替え
calc.beta0 <- function(res,mY=5,prob.beta=c(0.5,0.9),prob.delta=c(0.9,0.95),beta=1,delta=1,beta.est=TRUE,delta.est=FALSE,beta.range=c(0,1),delta.range=c(0.1,5),Fm2.max=5,thin=1,step1=0.2,tol=0.0001,
Btar=res$Btar, # いちおう、各種管理基準値は外からでも与えられるようにした
Blow=res$Blow,
Blim=res$Blim,
Bban=res$Bban,
Fmsy=res$Fmsy)
{
stockid <- res$stockid
farg <- res$farg
N <- res$N
GT <- res$GT
sim0 <- res$sim0
nyear <- res$nyear
lag <- res$lag
future.vpa1 <- get(res$future.function.name)
fm <- res$fm
w.recent <- res$w.recent
nY <- nyear
eyear <- res$eyear
SRdata <- res$SRdata
seed <- res$seed
B.cur <- SRdata$SSB[length(SRdata$SSB)]
# alpha & abc calculation
target.func <- res$target.func
targ <- res$TARres
Nlast <- targ$naa
error.last <- targ$eaa+w.recent
Dim.targ <- dim(Nlast)
if (lag == 0) Blast <- NULL else Blast <- targ$ssb[,Dim.targ[2]-(lag-1),]
targ.calc <- function(x) colSums(target.func(Fmsy,farg,mY=mY,eyear=1,seed=seed,N=N,naa0=Nlast[,Dim.targ[2],],eaa0=error.last[Dim.targ[2],],ssb0=Blast,beta=x,delta=0,Blim=Blim,Bban=Bban)$ssb[,1,])
if (beta.est){
beta.f1 <- function(x) {
ssb.tmp <- targ.calc(x)
Prob1 <- mean(ssb.tmp > Btar)
x1 <- Prob1-prob.beta[1]
dist1 <- x1^2
dist1
}
beta.f2 <- function(x) {
ssb.tmp <- targ.calc(x)
Prob2 <- mean(ssb.tmp > Blim)
x2 <- Prob2-prob.beta[2]
dist1 <- x2^2
dist1
}
res.beta1 <- optimize(beta.f1,beta.range)
res.beta2 <- optimize(beta.f2,beta.range)
beta <- min(res.beta1$minimum,res.beta2$minimum)
# beta <- floor(beta * 100)/100
}
## グラフによる図示
# ssb.msy <- apply(targ$ssb,c(2,3),sum)[5,]
# plot(density(ssb.msy),type="l",title="SSB")
# abline(v=Blim,lty=2,col=2)
# abline(v=Btar,lty=2)
# mean(ssb.msy>Blim)
# mean(ssb.msy>Btar)
ssb.tmp <- targ.calc(beta)
Prob.b1 <- mean(ssb.tmp > Btar)
Prob.b2 <- mean(ssb.tmp > Blim)
cat("beta=",round(beta,2),"; Prob SSB>Btar=",Prob.b1, "; Prob SSB>Blim=",Prob.b2,"; \n")
out <- list(beta=beta)
invisible(list(out,future.pred))
}
calc.beta2 <- function(MSY.input, # Fmsy(Ftarget)で漁獲するような将来予測の引数
mY=5,prob.beta=c(0.5,0.9),beta=1,delta=1,
beta.est=TRUE,beta.range=c(0,1),
Btar=res$Btar, # 各種管理基準値は外から与える
Blim=res$Blim,
Bban=res$Bban,
Fmsy=res$Fmsy)
{
stockid <- res$stockid
farg <- res$farg
N <- res$N
GT <- res$GT
sim0 <- res$sim0
nyear <- res$nyear
lag <- res$lag
future.vpa1 <- get(res$future.function.name)
fm <- res$fm
w.recent <- res$w.recent
nY <- nyear
eyear <- res$eyear
SRdata <- res$SRdata
seed <- res$seed
B.cur <- SRdata$SSB[length(SRdata$SSB)]
# alpha & abc calculation
target.func <- res$target.func
targ <- res$TARres
Nlast <- targ$naa
error.last <- targ$eaa+w.recent
Dim.targ <- dim(Nlast)
if (lag == 0) Blast <- NULL else Blast <- targ$ssb[,Dim.targ[2]-(lag-1),]
targ.calc <- function(x) colSums(target.func(Fmsy,farg,mY=mY,eyear=1,seed=seed,N=N,naa0=Nlast[,Dim.targ[2],],eaa0=error.last[Dim.targ[2],],ssb0=Blast,beta=x,delta=0,Blim=Blim,Bban=Bban)$ssb[,1,])
if (beta.est){
beta.f1 <- function(x) {
ssb.tmp <- targ.calc(x)
Prob1 <- mean(ssb.tmp > Btar)
x1 <- Prob1-prob.beta[1]
dist1 <- x1^2
dist1
}
beta.f2 <- function(x) {
ssb.tmp <- targ.calc(x)
Prob2 <- mean(ssb.tmp > Blim)
x2 <- Prob2-prob.beta[2]
dist1 <- x2^2
dist1
}
res.beta1 <- optimize(beta.f1,beta.range)
res.beta2 <- optimize(beta.f2,beta.range)
beta <- min(res.beta1$minimum,res.beta2$minimum)
# beta <- floor(beta * 100)/100
}
## グラフによる図示
# ssb.msy <- apply(targ$ssb,c(2,3),sum)[5,]
# plot(density(ssb.msy),type="l",title="SSB")
# abline(v=Blim,lty=2,col=2)
# abline(v=Btar,lty=2)
# mean(ssb.msy>Blim)
# mean(ssb.msy>Btar)
ssb.tmp <- targ.calc(beta)
Prob.b1 <- mean(ssb.tmp > Btar)
Prob.b2 <- mean(ssb.tmp > Blim)
cat("beta=",round(beta,2),"; Prob SSB>Btar=",Prob.b1, "; Prob SSB>Blim=",Prob.b2,"; \n")
out <- list(beta=beta)
invisible(list(out,future.pred))
}
############################# ここまで
### dynamics MSYを計算してみる
dyn.msy <- function(naa.past,naa.init=NULL,fmsy,a,b,resid,resid.year,waa,maa,M,SR=TRUE){
nyear <- length(resid)
if(is.null(naa.init)) nage <- nrow(naa.past) else nage <- length(naa.init)
naa <- matrix(0,nage,nyear)
ssb <- numeric()
if(is.null(naa.init)) naa[,1] <- naa.past[,colnames(naa.past)==min(resid.year)]
else naa[,1] <- naa.init
colnames(naa) <- resid.year
if(is.null(naa.init)){
waa <- waa[,colnames(naa.past)%in%resid.year]
maa <- maa[,colnames(naa.past)%in%resid.year]
M <- M[,colnames(naa.past)%in%resid.year]
}
for(i in 2:nyear){
ssb[i-1] <- sum(naa[,i-1]*waa[,i-1]*maa[,i-1],na.rm=T)
if(SR==TRUE){
naa[1,i] <- HS(ssb[i-1],a,b)*exp(resid[i])
}
else{
naa[1,i] <- naa.past[1,i]
}
for(j in 2:(nage-1)) naa[j,i] <- naa[j-1,i-1] * exp(-fmsy[j-1]-M[j-1,i-1])
naa[nage,i] <- naa[nage-1,i-1] * exp(-fmsy[j-1]-M[j-1,i-1]) + naa[nage,i-1] * exp(-fmsy[nage]-M[nage,i-1])
}
i <- nyear ; ssb[i] <- sum(naa[,i]*waa[,i]*maa[,i])
list(naa=naa,ssb=ssb)
}
plot.HCR <- function(beta=1,bban=0,blimit=1,btarget=2,add=FALSE,yscale=1.3,xlim=NULL,
Fmsy=1,scale=1,
ssb.cur=NULL,...) {
if(is.null(xlim)) xlim <- c(0,max(c(bban,balimit,btarget)))/scale
b.tmp <- seq(from=0,to=max(xlim),length=300)
y <- (b.tmp-bban)/(blimit-bban)*beta
y <- ifelse(b.tmp>blimit,beta,y)
y <- ifelse(y<0,0,y)
if(!isTRUE(add)) plot(b.tmp/scale,Fmsy*y,type="n",ylim=c(0,yscale),xlab="SSB",ylab="multiplier to current F",xlim=xlim/scale)
points(b.tmp/scale,Fmsy*y,type="l",...)
abline(h=Fmsy,col="gray")
text(0,Fmsy+0.02,paste("Fmsy=",round(Fmsy,2),"Fcurrent"),adj=0)
text(0,Fmsy*beta+0.02,paste("Ftarget=",round(beta*Fmsy,2),"Fcurrent (",round(beta,2),"*Fmsy)",sep=""),adj=0)
if(!is.null(ssb.cur)){
Frec <- (ssb.cur-bban)/(blimit-bban)
Frec <- ifelse(Frec<0,0,Frec)
Frec <- ifelse(Frec>1,1,Frec)
lines(c(0,ssb.cur/scale,ssb.cur/scale),c(Frec*beta*Fmsy,Frec*beta*Fmsy,0),lty=2)
# points(ssb.cur/scale,Frec*beta*Fmsy,lty=2,pch=4)
# text(0,0.8*beta*Fmsy+0.02,
# paste("F2018=",round(Frec*beta*Fmsy,2),"","Fcurrent (",
# round(Frec,2),"*",round(beta,2),"*",round(Fmsy,2),"*Fcurrent)",sep=""),adj=0)
}
}
draw.refline <- function(reftable,horiz=TRUE,scale=1000,lwd=3){
if(horiz==FALSE){
abline(v=reftable[1:4,2]/scale,
col=c("darkgreen","darkblue","darkred","red"),lwd=lwd,lty="22")
abline(v=reftable[1:4,1]/scale,
col=c("darkgreen","darkblue","darkred","red"),lwd=lwd,lty=1)
}
else{
abline(h=reftable[1:4,2]/scale,
col=c("darkgreen","darkblue","darkred","red"),lwd=lwd,lty="22")
abline(h=reftable[1:4,1]/scale,
col=c("darkgreen","darkblue","darkred","red"),lwd=lwd,lty=1)
}
}
#---- 期間内のRPSをサンプリング。平均値と中央値の差を補正するかどうか?など
RPS.simple.rec <- function(ssb,vpares,
rec.arg=list(rps.year=NULL, # 点推定値のrpsを計算する期間
upper.ssb=Inf, # 親魚資源量の上限(単位はトン?)
upper.recruit=Inf,
sample.year = NULL, # リサンプリング期間。rps.yearと異なる範囲を使う場合、設定する
bias.correction=TRUE, # stochasticのときに平均値と中央値の比率を使うもの)
rpsmean=FALSE),# 決定論的予測をRPSの平均でおこなうか、中央値でおこなうか。スケトウでは平均、その他の魚種は中央値。
deterministic=FALSE,rec.resample=NULL # ここは外から指定する必要ない
){
# argname <- c("rps.year","upper.ssb","upper.recruit","sample.year","bias.corrected","rpsmean")
# tmp <- !(names(rec.arg) %in% argname)
# if(sum(tmp)>0) stop(paste(names(rec.arg)[tmp]," no such arguments in RPS.simple.rec"))
if(is.null(rec.arg$bias.correction)) rec.arg$bias.correction <- TRUE
if(is.null(rec.arg$rpsmean)) rec.arg$rpsmean <- FALSE
if(is.null(rec.arg$rps.year)) rec.arg$rps.year <- as.numeric(colnames(vpares$naa))
if(is.null(rec.arg$sample.year)) rec.arg$sample.year <- rec.arg$rps.year
# browser()
names(rec.arg)
if(is.null(rec.resample)){
min.age <- min(as.numeric(rownames(vpares$ssb)))
if(min.age==0) slide.tmp <- TRUE else slide.tmp <- -1:-min.age
rps.data <- data.frame(year=years <- as.numeric(names(colSums(vpares$ssb,na.rm=T))),
ssb=SSB <- as.numeric(colSums(vpares$ssb,na.rm=T)),
recruit=rec <- as.numeric(c(vpares$naa[1,slide.tmp],
rep(NA,min.age))))
rps.data$rps <- rps <- rps.data$recruit/rps.data$ssb
rps.range <- as.numeric(rps[years %in% rec.arg$rps.year])
rps.med <- median(rps.range) # 点推定のためのrps
rps.mean <- mean(rps.range) # 点推定のためのrps
rec.resample <- as.numeric(rps[years %in% rec.arg$sample.year]) # リサンプリングのためのrps
# sample.mean <- mean(sample.range) # リサンプリングのためのrps
# sample.median <- median(sample.range) # リサンプリングのためのrps
# if(rec.arg$bias.correction==TRUE){
## rec.resample <- sample.range/rps.mean*rps.med
## rec.resample <- sample.range/sample.mean*rps.med # ここは本当にsample.meanで良いのか?(サンプル期間が同じ場合、sample.mean=rps=meanなので問題ない。期間が異なる場合、rps.meanを使うとsample.range/rps.meanの平均が1にならないため、やはりsample.meanを使うのが適切) => 分母はrps.medでなくsample.meanでは?
# rec.resample <- sample.range/sample.mean*sample.median
# }
# else{
# rec.resample <- sample.range
# }
}
ssb.tmp <- ifelse(ssb>rec.arg$upper.ssb,rec.arg$upper.ssb,ssb)
rec_determine <- ifelse(rec.arg$rpsmean,ssb.tmp * mean(rec.resample),ssb.tmp * median(rec.resample))
rec_random <- sample(rec.resample,length(ssb.tmp)-1,replace=TRUE) * ssb.tmp[-1]
if(rec.arg$bias.correction==TRUE){
rec_random <- rec_random * mean(rec.resample)/median(rec.resample)
}
rec <- c(rec_determine,rec_random)
rec2 <- ifelse(rec>rec.arg$upper.recruit,rec.arg$upper.recruit,rec)
# rec2 <- min(rec,rec.arg$upper.recruit)
return(list(rec=rec2,rec.resample=rec.resample))
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.