source("./tools/generate-testdata/rvpa1.9.4.r")
source("./tools/generate-testdata/future2.1.r")
#source("./tools/generate-testdata/future.r")
source("./tools/generate-testdata/utilities.r")
source("./tools/generate-testdata/stock_recruit.r") #stock_recruit.r 2019/11/26 ver.
caa <- read.csv(system.file("extdata","caa_pma.csv",package="frasyr"),row.names=1)
waa <- read.csv(system.file("extdata","waa_pma.csv",package="frasyr"),row.names=1)
maa <- read.csv(system.file("extdata","maa_pma.csv",package="frasyr"),row.names=1)
dat <- data.handler(caa=caa, waa=waa, maa=maa, M=0.5)
res_vpa_pma <- vpa(dat,fc.year=2009:2011,rec=585,rec.year=2011,tf.year = 2008:2010,
term.F="max",stat.tf="mean",Pope=TRUE,tune=FALSE,p.init=1.0)
SRdata_pma <- get.SRdata(res_vpa_pma)
SRmodel.list <- expand.grid(SR.rel = c("HS","BH","RI"), AR.type = c(0, 1), out.AR=c(TRUE,FALSE), L.type = c("L1", "L2"))
SR.list <- list()
for (i in 1:nrow(SRmodel.list)) {
SR.list[[i]] <- fit.SR(SRdata_pma, SR = SRmodel.list$SR.rel[i], method = SRmodel.list$L.type[i],
AR = SRmodel.list$AR.type[i], out.AR =SRmodel.list$out.AR[i], hessian = FALSE)
}
SRmodel.list$AICc <- sapply(SR.list, function(x) x$AICc)
SRmodel.list$delta.AIC <- SRmodel.list$AICc - min(SRmodel.list$AICc)
SR.list <- SR.list[order(SRmodel.list$AICc)] # AICの小さい順に並べたもの
(SRmodel.list <- SRmodel.list[order(SRmodel.list$AICc), ]) # 結果
SRmodel.base <- SR.list[[1]] # AIC最小モデルを今後使っていく
selectedSR <- sprintf("%s.recAR",SRmodel.base$input$SR[1])
# calc.yield curve=FALSEオプションの結果オブジェクトを作成(T/FでestMSYのtraceのオブジェクトサイズが違う)
# future fcurrent ----
res_future_Fcurrent_pma <- future.vpa(res_vpa_pma,
multi=1,
nyear=50, # 将来予測の年数
start.year=2012, # 将来予測の開始年
N=100, # 確率的計算の繰り返し回数
ABC.year=2013, # ABCを計算する年
waa.year=2009:2011, # 生物パラメータの参照年
maa.year=2009:2011,
M.year=2009:2011,
is.plot=TRUE, # 結果をプロットするかどうか
seed=1,
silent=TRUE,
recfunc=eval(parse(text=selectedSR))
, # 再生産関係の関数
# recfuncに対する引数
rec.arg=list(a=SRmodel.base$pars$a,b=SRmodel.base$pars$b,
rho=SRmodel.base$pars$rho, # ここではrho=0なので指定しなくてもOK
sd=SRmodel.base$pars$sd,resid=SRmodel.base$resid))
# est MSY----
res_MSY_pma <- est.MSY(res_vpa_pma, # VPAの計算結果
res_future_Fcurrent_pma$input, # 将来予測で使用した引数
seed=res_future_Fcurrent_pma$input$seed,
resid.year=0, # ARありの場合、最近何年分の残差を平均するかをここで指定する。ARありの設定を反映させたい場合必ずここを1以上とすること(とりあえず1としておいてください)。
N=100, # 確率的計算の繰り返し回数=>実際の計算では1000~5000回くらいやってください
calc.yieldcurve=FALSE,
PGY=c(0.95,0.9,0.6,0.1), # 計算したいPGYレベル。上限と下限の両方が計算される
onlylower.pgy=FALSE, # TRUEにするとPGYレベルの上限は計算しない(計算時間の節約になる)
B0percent=c(0.2,0.3,0.4),
Bempirical=c(round(tail(colSums(res_vpa_pma$ssb),n=1)),
round(max(colSums(res_vpa_pma$ssb))),
24000, # 現行Blimit
SRmodel.base$pars$b) # HSの折れ点
)
save(res_MSY_pma,file = "./inst/extdata/res_MSY_pma_calcyieldcurve_F.rda")
# 確率計算の繰り返し数を10のn乗として各々simN回計算
for(n in 2:4){
#n <-2
sample.size <- 10^n
# future fcurrent ----
res_future_Fcurrent_pma <- future.vpa(res_vpa_pma,
multi=1,
nyear=58, # 将来予測の年数
start.year=2012, # 将来予測の開始年
N=sample.size, # 確率的計算の繰り返し回数
ABC.year=2013, # ABCを計算する年
waa.year=2009:2011, # 生物パラメータの参照年
maa.year=2009:2011,
M.year=2009:2011,
is.plot=TRUE, # 結果をプロットするかどうか
seed=1,
silent=TRUE,
recfunc=eval(parse(text=selectedSR))
, # 再生産関係の関数
# recfuncに対する引数
rec.arg=list(a=SRmodel.base$pars$a,b=SRmodel.base$pars$b,
rho=SRmodel.base$pars$rho, # ここではrho=0なので指定しなくてもOK
sd=SRmodel.base$pars$sd,resid=SRmodel.base$resid))
# iteration num
simN <- 100
# est MSY----
for(i in 1:simN){
res_MSY_pma_test <- est.MSY(res_vpa_pma, # VPAの計算結果
res_future_Fcurrent_pma$input, # 将来予測で使用した引数
#seed=res_future_Fcurrent_pma$input$seed,
seed=i,
resid.year=0, # ARありの場合、最近何年分の残差を平均するかをここで指定する。ARありの設定を反映させたい場合必ずここを1以上とすること(とりあえず1としておいてください)。
N=sample.size, # 確率的計算の繰り返し回数=>実際の計算では1000~5000回くらいやってください
calc.yieldcurve=FALSE,
PGY=c(0.95,0.9,0.6,0.1), # 計算したいPGYレベル。上限と下限の両方が計算される
onlylower.pgy=FALSE, # TRUEにするとPGYレベルの上限は計算しない(計算時間の節約になる)
B0percent=c(0.2,0.3,0.4),
Bempirical=c(round(tail(colSums(res_vpa_pma$ssb),n=1)),
round(max(colSums(res_vpa_pma$ssb))),
24000, # 現行Blimit
SRmodel.base$pars$b) # HSの折れ点
)
#outfile <- sprintf("./tools/generate-testdata/res_MSY_pma_%d.rda",i)
#save(res_MSY_pma,file = outfile)
assign(sprintf("res_MSY_pma%d",i),res_MSY_pma_test)
}
# 結果meanとsdのオブジェクトサイズ確保
load(system.file("extdata","res_MSY_pma_calcyieldcurve_F.rda",package = "frasyr"))
res_MSY_pma_mean <- res_MSY_pma
res_MSY_pma_mean$summary <- c()
res_MSY_pma_mean$summaryAR <- c()
res_MSY_pma_mean$all.stat <- c()
res_MSY_pma_mean$all.statAR <- c()
res_MSY_pma_sd <- res_MSY_pma
res_MSY_pma_sd$summary <- c()
res_MSY_pma_sd$summaryAR <- c()
res_MSY_pma_sd$all.stat <- c()
res_MSY_pma_sd$all.statAR <- c()
#オブジェクト初期化
for(j in 1:length(names(res_MSY_pma$summary_tb))){
if(is.character(res_MSY_pma$summary_tb[[j]])){
res_MSY_pma_mean$summary_tb[[j]] <- res_MSY_pma_sd$summary_tb[[j]] <- res_MSY_pma$summary_tb[[j]]
}else if(is.logical(res_MSY_pma$summary_tb[[j]])){
res_MSY_pma_mean$summary_tb[[j]] <- res_MSY_pma_sd$summary_tb[[j]] <- res_MSY_pma$summary_tb[[j]]
}else{
res_MSY_pma_mean$summary_tb[[j]] <- res_MSY_pma_sd$summary_tb[[j]] <- res_MSY_pma_mean$summary_tb[[j]]*0
}
}
for(j in 1:length(res_MSY_pma$F.msy)){
res_MSY_pma_mean$F.msy[[j]] <- res_MSY_pma_sd$F.msy[[j]] <- 0
}
for(j in 1:length(names(res_MSY_pma$all.stat_tb))){
if(is.character(res_MSY_pma$all.stat_tb[[j]])){
res_MSY_pma_mean$all.stat_tb[[j]] <- res_MSY_pma_sd$all.stat_tb[[j]] <- res_MSY_pma$all.stat_tb[[j]]
}else if(is.logical(res_MSY_pma$all.stat_tb[[j]])){
res_MSY_pma_mean$all.stat_tb[[j]] <- res_MSY_pma_sd$all.stat_tb[[j]] <- res_MSY_pma$all.stat_tb[[j]]
}else{
res_MSY_pma_mean$all.stat_tb[[j]] <- res_MSY_pma_sd$all.stat_tb[[j]] <- res_MSY_pma_mean$all.stat_tb[[j]]*0
}
}
for(j in 1:length(names(res_MSY_pma$trace))){
res_MSY_pma_mean$trace[[j]] <- res_MSY_pma_sd$trace[[j]] <- c(0)
}
for(j in 1:length(res_MSY_pma$ssb.ar.mean)){
res_MSY_pma_mean$ssb.ar.mean[j] <- res_MSY_pma_sd$ssb.ar.mean[j] <- c(0)
}
res_MSY_pma_mean$SPR.msy <- res_MSY_pma_sd$SPR.msy <- 0
#for(i in 1:simN){
# readlile <- sprintf("./tools/generate-testata/res_MSY_pma%d.rda",i)
# assign(sprintf("res_vpa_pma%d",i), res_vpa_pma)
# }
checkitem <-c("summary_tb","all.stat_tb","trace")
for(l in 1:length(checkitem)){
tmp <- list()
for(i in 1:simN){
tmp[[i]] <- as.data.frame(eval(parse(text=paste(sprintf("res_MSY_pma%d$%s",i,checkitem[l])))))
}
for(k in 1:length(eval(parse(text=paste(sprintf("res_MSY_pma$%s",checkitem[l])))))){
each.tmp <- c()
for(i in 1:simN){
each.tmp <- cbind(each.tmp,tmp[[i]][k][,])
}
#print(sprintf("res_MSY_pma$%s",checkitem[l]))
if(is.numeric(tmp[[i]][k][,])){
mean.tmp <- c()
sd.tmp <- c()
for(j in 1:length(tmp[[i]][k][,])){
mean.tmp <- rbind(mean.tmp,mean(each.tmp[j,]))
sd.tmp <- rbind(sd.tmp, sd(each.tmp[j,]))
}
if(l == 1){
res_MSY_pma_mean$summary_tb[k] <- mean.tmp
res_MSY_pma_sd$summary_tb[k] <- sd.tmp
} else if(l == 2){
res_MSY_pma_mean$all.stat_tb[k] <- mean.tmp
res_MSY_pma_sd$all.stat_tb[k] <- sd.tmp
} else if(l == 3){
res_MSY_pma_mean$trace[k] <- mean.tmp
res_MSY_pma_sd$trace[k] <- sd.tmp
}
}
}
}
checkitem <-c("F.msy")
for(l in 1:length(checkitem)){
tmp <- list()
each.tmp <-c()
for(i in 1:simN){
tmp[[i]] <- as.data.frame(eval(parse(text=paste(sprintf("res_MSY_pma%d$%s",i,checkitem[l])))))
if(is.numeric(tmp[[i]][,])){
each.tmp <- cbind(each.tmp,tmp[[i]][,])
}
}
mean.tmp <- c()
sd.tmp <- c()
for(j in 1:length(tmp[[i]][,])){
mean.tmp <- rbind(mean.tmp,mean(each.tmp[j,]))
sd.tmp <- rbind(sd.tmp, sd(each.tmp[j,]))
}
res_MSY_pma_mean$F.msy <- mean.tmp
res_MSY_pma_sd$F.msy <- sd.tmp
}
checkitem <-c("ssb.ar.mean")
for(l in 1:length(checkitem)){
tmp <- list()
for(i in 1:simN){
tmp[[i]] <- as.data.frame(eval(parse(text=paste(sprintf("res_MSY_pma%d$%s",i,checkitem[l])))))
}
mean.mat <- c()
sd.mat <- c()
for(k in 1:ncol(eval(parse(text=paste(sprintf("res_MSY_pma$%s",checkitem[l])))))){
each.tmp <-c()
for(i in 1:simN){
each.tmp <- cbind(each.tmp,tmp[[i]][k][,])
}
mean.tmp <- c()
sd.tmp <- c()
for(j in 1:nrow(tmp[[i]][k])){
mean.tmp <- rbind(mean.tmp,mean(each.tmp[j,]))
sd.tmp <- rbind(sd.tmp, sd(each.tmp[j,]))
}
mean.mat <- cbind(mean.mat,mean.tmp)
sd.mat <- cbind(sd.mat, sd.tmp)
}
res_MSY_pma_mean$ssb.ar.mean <- mean.mat
res_MSY_pma_sd$ssb.ar.mean <- sd.mat
}
checkitem <-c("SPR.msy")
tmp <- list()
for(i in 1:simN){
tmp[[i]] <- as.data.frame(eval(parse(text=paste(sprintf("res_MSY_pma%d$%s",i,checkitem)))))
}
each.tmp <-c()
for(i in 1:simN){
each.tmp <- cbind(each.tmp,tmp[[i]][,])
}
mean.tmp <- mean(each.tmp)
sd.tmp <- sd(each.tmp)
res_MSY_pma_mean$SPR.msy <-mean.tmp
res_MSY_pma_sd$SPR.msy <-sd.tmp
outmeanfile <- sprintf("./inst/extdata/res_MSY_pma_mean_samplesize%d.rda",10^n)
outsdfile <- sprintf("./inst/extdata/res_MSY_pma_sd_samplesize%d.rda",10^n)
save(res_MSY_pma_mean,file=outmeanfile)
save(res_MSY_pma_sd,file=outsdfile)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.