R/stock_recruit.r

Defines functions hmm_SR fit.SR_pen fit.SR_tol tryall_SR check_consistent_w boot_steepness calc_steepness corSR bootSR.ggplot specify.Mesnil.gamma check.SRfit out.SR prof.likSR jackknife.SR bootSR.plot autocor.plot calc.residAR check.SRdist calc.StdResid fit.SRregime prof.lik boot.SR fit.SR2 fit.SRalpha fit.SR validate_sr get.SRdata

Documented in autocor.plot boot.SR bootSR.ggplot bootSR.plot boot_steepness calc.residAR calc.StdResid calc_steepness check_consistent_w check.SRdist check.SRfit corSR fit.SR fit.SR_pen fit.SRregime fit.SR_tol get.SRdata hmm_SR jackknife.SR out.SR prof.lik prof.likSR specify.Mesnil.gamma tryall_SR validate_sr

#'
#' @import magrittr
#' @import dplyr
#' @import tidyr
#' @import tibble
#' @importFrom magrittr %>%
#' @importFrom dplyr filter
#' @importFrom dplyr mutate
#' @importFrom dplyr group_by
#' @importFrom dplyr select

#'
NULL

#' VPAの結果から再生産関係推定用のデータを作成する
#'
#' VPA結果に放流尾数が含まれている場合、天然加入尾数=加入尾数ー放流尾数として計算する
#'
#' @param vpares VPAの結果のオブジェクト
#' @param R.dat VPAの結果オブジェクトでなく直接データを与えたい場合の加入の値
#' @param SSB.dat VPAの結果オブジェクトでなく直接データを与えたい場合の親魚量の値
#' @param release.dat VPAの結果オブジェクトでなく直接データを与えたい場合の放流尾数の値
#' @param return.df データフレームとして結果を返すか。このオプション関係なくデータフレームとして返すように変更したので、近いうちに廃止予定
#' @param weight.year fit.SRに渡すとき、フィットの対象とする年を指定する。特例として0を与えると、全部の年のデータが使われるようになる。返り値にweightという列が加わる。
#'
#' @encoding UTF-8
#'
#' @export
get.SRdata <- function(vpares=NULL,
                       R.dat = NULL,
                       SSB.dat = NULL,
                       release.dat=NULL,
                       years = NULL,
                       weight.year = NULL,
                       weight.data = NULL,
                       return.df = TRUE){

    is.release <- !is.null(release.dat) | !is.null(vpares$input$dat$release.alive)
    if(is.null(years) && !is.null(vpares)) years <- as.numeric(colnames(vpares$naa))

    # R.datとSSB.datだけが与えられた場合、それを使ってシンプルにフィットする
    if (!is.null(R.dat) & !is.null(SSB.dat)) {
        assertthat::assert_that(length(R.dat)==length(SSB.dat))
        if(!is.null(release.dat)){
            assertthat::assert_that(length(R.dat)==length(release.dat))
            R.dat <- R.dat - release.dat
        }
        if(is.null(years)){
            years <- 1:length(R.dat)
        }
        else{
            assertthat::assert_that(length(R.dat)==length(years))
        }
        dat <- tibble(R = R.dat,SSB = SSB.dat,year = years)
        if(!is.null(release.dat)) R.dat$release <- release.dat
    }

    # vpa dataが与えられた場合にはデータの整形をする
    if(!is.null(vpares)) {
        n <- ncol(vpares$naa)
        L <- as.numeric(rownames(vpares$naa)[1])

        dat      <- list()
        dat$SSB  <- as.numeric(colSums(vpares$ssb,na.rm = TRUE))
        dat$year <- as.numeric(colnames(vpares$ssb))
        dat$R   <-  as.numeric(vpares$naa[1,])

        # 加入年齢分だけずらす
        dat$R    <- dat$R[(L+1):n]
        dat$SSB  <- dat$SSB[1:(n-L)]
        dat$year <- dat$year[(L+1):n]

        dat <- as_tibble(dat)

        if(!is.null(vpares$input$dat$release.dat)) vpares$input$dat$release.alive <- vpares$input$dat$release.dat
        if(!is.null(vpares$input$dat$release.alive)){

          tmpfunc <- function(rdat, data_name){
            if(!is.null(rdat)){
              dat <- rdat[rownames(rdat) %in% as.character(L), ]
              rownames(dat) <- data_name
              dat %>% t() %>% as.data.frame() %>% rownames_to_column(var="year") %>%
                mutate(year=as.numeric(year)) %>% as_tibble()
            }
          }

          dat <- left_join(dat, tmpfunc(vpares$input$dat$release.alive, "release_alive"))
          dat <- dat %>% mutate(R=R-release_alive)
          if(!is.null(vpares$input$dat$release.all)){
            dat <- left_join(dat, tmpfunc(vpares$input$dat$release.all, "release_all"))
            dat <- left_join(dat, tmpfunc(vpares$input$dat$release.ratealive, "release_ratealive"))
          }
        }

        dat$R[dat$R<0] <- 0.001

        # データの抽出
        dat <- dat %>% dplyr::filter(year%in%years)
    }

    assertthat::assert_that(all(dat[["R"]] > 0))

    #if (return.df == TRUE) return(data.frame(year = dat$year,
    #SSB  = dat$SSB,
    #                                         R    = dat$R,
    #                                             release=dat$release))
    #return(dat[c("year","SSB","R","release")])

#    dat.df <- data.frame(year = dat$year,
#                         SSB  = dat$SSB,
#                         R    = dat$R)
#    if(is.release){
#        dat.df$release_alive <- dat$release_alive
#        if(!is.null(dat$release_all)) dat.df$release_all <- dat$release_all
#        if(!is.null(dat$release_ratealive))  dat.df$release_ratealive <- dat$
#    }

  if(!is.null(weight.year)){
    dat <- dat %>% mutate(weight=0)
    dat$weight[dat$year %in% weight.year]  <- 1

    if(length(weight.year)==1 && weight.year==0){
      dat <- dat %>% mutate(weight=1)
    }
  }
  if(is.null(weight.year)){
    dat <- dat %>% mutate(weight=1)
  }

  if(!is.null(weight.data)){
    assertthat::assert_that(all(dat$year==weight.data$year))
    if("weight" %in% names(dat)) dat <- dat %>% select(-weight)
    dat <- dat %>% left_join(weight.data)
  }

  dat <- dat %>% select(year, SSB, R, weight, everything())
  return(dat)

}

#' 再生産関係のパラメタが想定した形になっているかを確認
#'
#' @inheritParams fit.SR
#' @param res_sr fit.SR の結果
validate_sr <- function(SR = NULL, method = NULL, AR = NULL, out.AR = NULL, res_sr = NULL) {
  if (!is.null(SR)) {
    assertthat::assert_that(
      length(SR) == 1,
      SR %in% c("HS", "BH", "RI","Mesnil")
    )
  }
  if (!is.null(method)) {
    method = as.character(method)
    assertthat::assert_that(method %in% c("L1", "L2"))
  }
  if (!is.null(AR)) {
    assertthat::assert_that(
      is.numeric(AR),
      AR %in% c(0, 1)
    )
  }
  if (!is.null(out.AR)) {
    assertthat::assert_that(is.logical(out.AR))
  }
  if (!is.null(res_sr)) {
    validate_sr(SR     = res_sr$input$SR,
                method = res_sr$input$method,
                AR     = res_sr$input$AR,
                out.AR = res_sr$input$out.AR)
  }
}

#' 再生産関係の推定
#'
#' 3種類の再生産関係の推定を、最小二乗法か最小絶対値法で、さらに加入の残差の自己相関を考慮して行うことができる
#' @param SRdata \code{get.SRdata}で作成した再生産関係データ
#' @param SR 再生産関係 (\code{"HS"}: Hockey-stick, \code{"BH"}: Beverton-Holt, \code{"RI"}: Ricker, \code{"Mesnil"}: Continuous HS)
#' @param method 最適化法(\code{"L2"}: 最小二乗法, \code{"L1"}: 最小絶対値法)
#' @param AR 自己相関を推定するか(1), しないか(0)
#' @param out.AR 自己相関係数を一度再生産関係を推定したのちに、外部から推定するか(1), 内部で推定するか(0)
#' @param length 初期値を決める際のgridの長さ
#' @param p0 \code{optim}で設定する初期値
#' @param bio_par data.frame(waa=c(100,200,300),maa=c(0,1,1),M=c(0.3,0.3,0.3)) のような生物パラメータをあらわすデータフレーム。waaは資源量を計算するときのweight at age, maaはmaturity at age, Mは自然死亡係数。これを与えると、steepnessやR0も計算して返す
#' @param plus_group hなどを計算するときに、プラスグループを考慮するか
#' @param HS_fix_b SR=HSのとき、折れ点bを固定して計算したい時に値を代入。デフォルットはNULL。ただし、min(SSB)/100より小さい値は入れないこと。
#' @param gamma \code{SR="Mesnil"}のときに使用するsmoothing parameter
#'
#' @encoding UTF-8
#' @examples
#' \dontrun{
#' data(res_vpa)
#' SRdata <- get.SRdata(res_vpa)
#' resSR <- fit.SR(SRdata, SR = c("HS","BH","RI")[1],
#'                 method = c("L1","L2")[2], AR = 1,
#'                 out.AR = TRUE)
#' resSR$pars
#'
#' # When outputting steepness, create a bio_par object using derive_biopar function with the res_vpa object and the corresponding year as its argument, and put the bio_par object in the argument of fit.SR.
#' bio_par <- derive_biopar(res_obj=res_vpa,derive_year = 2010)
#' resSR <- fit.SR(SRdata, SR = c("HS","BH","RI")[1],
#'                 method = c("L1","L2")[2], AR = 1,
#'                 out.AR = TRUE,bio_par=bio_par)
#' resSR$pars
#' resSR$steepness
#' }
#' @return 以下の要素からなるリスト
#' \describe{
#' \item{\code{input}}{使用した引数のリスト}
#' \item{\code{pars}}{推定されたパラメータ}
#' \item{\code{opt}}{\code{optim}の結果オブジェクト}
#' \item{\code{resid}}{再生産関係から予測値からの加入量の残差}
#' \item{\code{resid2}}{自己相関のを推定したうえでの加入の残差(自己相関なしの時\code{resid}と等しくなる)}
#' \item{\code{loglik}}{対数尤度}
#' \item{\code{k}}{推定したパラメータ数}
#' \item{\code{AIC}}{AIC (\code{out.AR=TRUE}のときは自己相関推定前の結果)}
#' \item{\code{AICc}}{AICc (\code{out.AR=TRUE}のときは自己相関推定前の結果)}
#' \item{\code{BIC}}{BIC (\code{out.AR=TRUE}のときは自己相関推定前の結果)}
#' \item{\code{AIC.ar}}{\code{out.AR=TRUE}のときに\code{acf}関数で得られた自己相関を推定しない場合(0)とする場合(1)のAICの差}
#' \item{\code{pred}}{予測された再生産関係}
#' \item{\code{steepness}}{bio_parを与えたときに、steepness (h) やR0(漁獲がない場合の平均加入尾数), SB0(漁獲がない場合の平均親魚量)なども追加的に返す}
#' }
#'
#' @export
#'

fit.SR <- function(SRdata,
                   SR="HS",
                   method="L2",
                   AR=1,
                   # TMB=FALSE,
                   hessian=FALSE,w=NULL,
                   length=20,
                   max.ssb.pred=1.3, # 予測値を計算するSSBの最大値(観測された最大値への乗数)
                   p0=NULL,
                   out.AR = TRUE, #自己相関係数rhoを外で推定するか
                   bio_par = NULL,
                   plus_group = TRUE,
                   is_jitter = FALSE,
                   HS_fix_b = NULL,
                   gamma=0.01
){
  validate_sr(SR = SR, method = method, AR = AR, out.AR = out.AR)

  argname <- ls()
  arglist <- lapply(argname,function(xx) eval(parse(text=xx)))
  names(arglist) <- argname

  tmp <- check_consistent_w(w, SRdata)
  SRdata <- tmp$SRdata
  w <- arglist$w <- tmp$w


  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)
  if (SR=="Mesnil") SRF <- function(x,a,b) 0.5*a*(x+sqrt(b^2+gamma^2/4)-sqrt((x-b)^2+gamma^2/4))

  if (length(SRdata$R) != length(w)) stop("The length of 'w' is not appropriate!")

  if (!is.null(HS_fix_b) && SR!="HS" ) stop("Set SR type as HS if you use HS_fix_b option!")
  if (!is.null(HS_fix_b) && HS_fix_b < (min(ssb)/100) ) stop("Set larger HS_fix_b than min(SSB)/100 if you use HS_fix_b option!")

  one_max = max(SRdata$year[w>0])
  zero_min =ifelse(sum(w==0)>0, min(SRdata$year[w==0]),one_max)

  if (method == "L2" && AR==1 && out.AR==FALSE && zero_min<one_max) { # For Jackknife

    obj.f = function(a,b,rho,out = "nll") {
      w2 = as.numeric(w>0)
      before_zero = rep(0,N)
      for (i in 1:N) {
        if (w2[i]>0) {
          if (i == 1) {
            w2[i] <- w[i]*(1-rho^2)
          } else {
            for (j in 1:(i-1)) {
              if (rev(w[1:(i-1)])[j] == 0) {
                before_zero[i] = before_zero[i]+1
                w2[i] <- w2[i]+rho^(2*j)
              } else break
            }
            if (before_zero[i] == i-1) {
              w2[i] <- w[i]*(1-rho^2)
            } else {
              w2[i] <- w[i]/w2[i]
            }
          }
        }
      }

      resid <- sapply(1:N,function(i) log(rec[i]) - log(SRF(ssb[i],a,b)))
      pred_resid <- NULL
      for (i in 1:N) {
        if (i==1 || before_zero[i] == i-1) {
          pred_resid <- c(pred_resid,0)
        } else {
          pred_resid <- c(pred_resid,rho^(1+before_zero[i])*resid[i-1-before_zero[i]])
        }
      }
      sd2 = sum(w2*(resid-pred_resid)^2)/sum(w2) #variance
      sd = sqrt(sd2) #SD

      sigma = NULL
      for (i in 1:N) {
        if (w2[i]==0) {
          sigma = c(sigma,1000)
        } else {
          sigma = c(sigma,sqrt(sd2/w2[i]))
        }
      }

      nll <- -sum(w*dnorm(resid,pred_resid,sigma,log=TRUE))
      if (out=="nll") return(nll)
      if (out=="resid") return(resid)
      if (out=="resid2") return(resid-pred_resid)
      if (out=="sd") return(sd)
    }
  } else {
    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(HS_fix_b)){
      if (is.null(p0)) {
    a.range <- range(rec/ssb)
    b.range <- range(1/ssb)
    if (SR == "HS" | SR=="Mesnil") 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" | SR =="Mesnil",-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" | SR == "Mesnil") {
      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)
    #if (rep.opt) {
    for (i in 1:100) {
      if(is_jitter==FALSE) par2 <- opt$par
      if(is_jitter==TRUE)  par2 <- opt$par + rnorm(length(opt$par),0,5)
      opt2 <- optim(par2,obj.f2)
      if (abs(opt$value-opt2$value)<1e-6) break
      opt <- opt2
    }
    #}
    opt <- optim(opt$par,obj.f2,method="BFGS",hessian=hessian)

  }
  else{
    if (is.null(p0)) {
      a.range <- c(min(rec/ssb),max(rec)/(min(ssb)/100))
      grids <- as.matrix(seq(a.range[1],a.range[2],len=length))
      init <- as.numeric(grids[which.min(sapply(1:length(grids),function(i) obj.f(grids[i,1],HS_fix_b,0))),])
      init[1] <- log(init[1])
      if (AR != 0 && !isTRUE(out.AR)) init[2] <- 0
    } else {
      init = p0
    }

    if (AR == 0 || out.AR) {
      obj.f2 <- function(x) obj.f(exp(x[1]),HS_fix_b,0)
      opt <- optim(init,obj.f2,method ="Brent",upper = log(max(a.range)), lower = log(min(a.range)))
      #if (rep.opt) {
      #}

    } else {
      obj.f2 <-  function(x) obj.f(exp(x[1]),HS_fix_b,1/(1+exp(-x[2])))

      opt <- optim(init,obj.f2)
      #if (rep.opt) {
      for (i in 1:100) {
        if(is_jitter==FALSE) par2 <- opt$par
        if(is_jitter==TRUE)  par2 <- opt$par + rnorm(length(opt$par),0,5)
        opt2 <- optim(par2,obj.f2)
        if (abs(opt$value-opt2$value)<1e-6) break
        opt <- opt2
      }
      #}
      opt <- optim(opt$par,obj.f2,method="BFGS",hessian=hessian)

    }

  }

  Res <- list()
  Res$input <- arglist
  Res$obj.f <- obj.f
  Res$obj.f2 <- obj.f2
  Res$opt <- opt

  if(is.null(HS_fix_b)){
    a <- exp(opt$par[1])
    b <- ifelse(SR=="HS"|SR=="Mesnil",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]))))
  } else{
    a <- exp(opt$par[1])
    b <- HS_fix_b
    rho <- ifelse(AR==0,0,ifelse(out.AR,0,1/(1+exp(-opt$par[2]))))
  }

  if (method == "L2" && AR==1 && out.AR==FALSE && zero_min<one_max) {
    resid = obj.f(a=a,b=b,rho=rho,out="resid")
    resid2 = obj.f(a=a,b=b,rho=rho,out="resid2")
    sd <- sd.pred <- obj.f(a=a,b=b,rho=rho,out="sd")
  } else {
    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 <- sd.pred <- sqrt(rss/NN)
    # } else {
      if (method=="L1") {
        rss <- w[1]*abs(resid2[1])*sqrt(1-rho^2)
        for(i in 2:N) rss <- rss + w[i]*abs(resid2[i])
        sd.pred <- sum(abs(w*resid2))/NN
        sd.pred <- sqrt(2)*sd.pred
    }
    # 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$sd.pred = sd.pred
  # if (sd.obs) {
    # if (AR==0 && method != "L2") sd = sqrt(sum(resid[w==1]^2)/sum(w))
  # }
  Res$pars <- c(a,b,sd,rho)

  if (method!="L2") {
    if (AR!=0) {
      if (!isTRUE(out.AR)) {
        message("L1 & out.AR=FALSE is NOT recommended")
      }
      arres <- ar(resid,aic=FALSE,order.max=1,demean=FALSE,method="mle")
      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] <- Res$sd.pred <-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)

  if(!is.null(bio_par)){
    if(SR!="Mesnil") Res$steepness <- calc_steepness(SR=SR,Res$pars,bio_par$M,bio_par$waa,bio_par$maa,plus_group=plus_group) else{ #add gamma to Res$pars if SR=Mesnil
      pars_and_gamma <- data.frame(Res$pars,gamma)
      Res$steepness <- calc_steepness(SR=SR,pars_and_gamma,bio_par$M,bio_par$waa,bio_par$maa,plus_group=plus_group)
    }
  }

  class(Res) <- "fit.SR"

  return(Res)
}



#' L1とL2のmixtureによる再生産関係の推定(beta版なのでexportしない)
#'
#' 3種類の再生産関係の推定を、最小二乗法か最小絶対値法で、さらに加入の残差の自己相関を考慮して行うことができる
#' @param SRdata \code{get.SRdata}で作成した再生産関係データ
#' @param SR 再生産関係 (\code{"HS"}: Hockey-stick, \code{"BH"}: Beverton-Holt, \code{"RI"}: Ricker)
#' @param alpha alpha:(1-alpha)の比でL1:L2を混ぜる
#' @param length 初期値を決める際のgridの長さ
#' @param rep.opt \code{TRUE}で\code{optim}による最適化を収束するまで繰り返す
#' @param p0 \code{optim}で設定する初期値
#' @encoding UTF-8
# #' @export
#' @noRd
#'
fit.SRalpha <- function(SRdata,
                   SR="HS",
                   alpha=0,
                   hessian=FALSE,
                   w=rep(1,length(SRdata$R)),
                   length=20,
                   max.ssb.pred=1.3, # 予測値を計算するSSBの最大値(観測された最大値への乗数)
                   p0=NULL,
                   rep.opt = TRUE
){

  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を計算するサンプル数

  validate_sr(SR = SR)
  #  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,out="nll"){ #rhoは無し
    resid <- sapply(1:N,function(i) log(rec[i]) - log(SRF(ssb[i],a,b)))

    # L2 part
    rss_L2 = sum(w*resid^2)
    sd = sqrt(rss_L2/NN)
    obj_L2 = -sum(w*dnorm(resid,0,sd,log=TRUE))

    # L1 part
    rsa_L1 = sum(w*abs(resid))
    phi = rsa_L1/NN
    obj_L1 = -sum(w*(-log(2*phi)-abs(resid/phi)))

    obj = alpha*obj_L2 + (1-alpha)*obj_L1
    SD = sqrt(alpha*sd^2 + (1-alpha)*2*phi^2)
    if (out=="nll") return(obj)
    if (out=="sd") return(SD)
    if (out=="resid") return(resid)
  }

  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]))),])
    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") {
    obj.f2 <- function(x) obj.f(exp(x[1]),min(ssb)+(max(ssb)-min(ssb))/(1+exp(-x[2])))
  } else {
    obj.f2 <- function(x) obj.f(exp(x[1]),exp(x[2]))
  }

  opt <- optim(init,obj.f2)
  if (rep.opt) {
    for (i in 1:100) {
      opt2 <- optim(opt$par,obj.f2)
      if (abs(opt$value-opt2$value)<1e-6) break
      opt <- opt2
    }
  }
  opt <- optim(opt$par,obj.f2,method="BFGS",hessian=hessian)

  Res <- list()
  Res$input <- arglist
  Res$obj.f <- obj.f
  Res$obj.f2 <- obj.f2
  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]))))
  rho <- 0

  sd = obj.f(a=a,b=b,out="sd")
  resid = obj.f(a=a,b=b,out="resid")

  Res$resid <- resid
  Res$resid2 <- resid
  Res$pars <- c(a,b,sd,rho)

  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)

  class(Res) <- "fit.SRalpha"
  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

  validate_sr(SR = SR, method = method, AR = AR)
  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)
}

#' 再生産関係のブートストラップ
#'
#' ①残差のパラメトリックブートストラップ、②残差のノンパラメトリックブートストラップ(リサンプリング)、③データのブートストラップ(リサンプリング)が行える
#'
#' 残差のパラメトリックブートストラップはmethod="p"で実行可能で、自己相関なしの場合、fit.SRのpars$sdを標準偏差とする正規分布からランダムに乱数を発生させ、予測値からのずれを加えて加入量のブートストラップデータを生成し、再推定している
#'
#' 自己相関ありの場合、fit.SRのpars$sdを1から自己相関係数rho^2を引き2乗根をとって除したもの(pars$sd/sqrt(1-pars$rho^2))を標準偏差とする正規分布からランダムに乱数(epsilon_t)を発生させ、 毎年の残差(resid_t) を resid_(t+1) = rho × resid_t + epsilon_t とし、resid_tを予測値に加えたブートストラップデータを生成し、再推定している
#'
#' 自己相関ありの場合はノンパラメトリックブートストラップは使わずにパラメトリックブートストラップを用いること
#'
#' 残差のノンパラメトリックブートストラップはmethod="n"で実行可能で、残差の確率分布を仮定せず、残差を重複ありでリサンプリングして、加入量のブートストラップデータを生成する
#'
#' データのブートストラップはmethod="d"で実行可能で、データを重複ありでリサンプリングしたデータを使用して、再生産関係の再推定を行う
#'
#' 親魚量データもリサンプリングにより変化するため、親魚量の不確実性も考慮されることになることから、親魚量データに偏りがあったり、データ数が少なかったり、あるデータ点に推定値が大きく依存している場合はバイアスや不確実性が大きくなりやすい
#'
#' 図のプロットにはbootSR.plotを使用する
#'
#' 自己相関を推定していない場合は最後のrhoの図は表示されない
#'
#' fit.SRの引数にbio_parを入れてスティープネスを計算した場合、SB0、R0、B0、hの図もされる
#'
#' 図の出力bootSR.plotのオプションoutput=Tとすると、各パラメータのヒストグラムが出力される
#'
#' @import purrr
#' @param Res \code{fit.SR}か\code{fit.SRregime}のオブジェクト
#' @param method 残差ブートストラップ(パラメトリック ("p") かノンパラメトリック ("n")) もしくはデータブートストラップ("d")
#' @param n ブートストラップの回数(例では100回だが500回あるいは1000回を推奨)
#' @encoding UTF-8
#' @examples
#' \dontrun{
#' data(res_vpa)
#' SRdata <- get.SRdata(res_vpa)
#' bio_par <- derive_biopar(res_obj=res_vpa,derive_year = 2010)
#' resL1outer = fit.SR(SRdata, SR = "HS", method = "L1", out.AR = TRUE, AR = 1,bio_par=bio_par)
#'
#' # example if parametric bootstrap
#' boot.res1 = boot.SR(resL1outer, n = 100, method = "p")
#' bootSR.plot(boot.res1)
#'
#' # example if non-parametric bootstrap
#' boot.res2 = boot.SR(resL1outer, n = 100, method = "n")
#' bootSR.plot(boot.res2)
#'
#' #' # example if data bootstrap (optional method)
#' boot.res3 = boot.SR(resL1outer, n = 100, method = "d")
#' bootSR.plot(boot.res3)

#' }
#' @export
#'
boot.SR <- function(Res,method="p",n=100,seed=1){

  argname <- ls()
  arglist <- lapply(argname,function(xx) eval(parse(text=xx)))
  names(arglist) <- argname
  N <- length(Res$input$SRdata$SSB)

  validate_sr(res_sr = Res)
  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)

  set.seed(seed)
  if (class(Res) == "fit.SR") { #fit.SR
    N <- length(Res$input$SRdata$SSB)
    RES = lapply(1:n, function(j){
        sd <- sapply(1:N, function(i) ifelse(i==1,Res$pars$sd/sqrt(1-Res$pars$rho^2),Res$pars$sd))
        for (j in 1:100) {
          if (method=="d") { #data bootstrap
            Year0 = Res$input$SRdata$year[Res$input$w==1]
            Year_r = sample(Year0,replace=TRUE)
            w_r = purrr::map_dbl(Year0, function(x) sum(Year_r==x))
            res.b <- Res
            res.b$input$w <- res.b$input$SRdata$weight <-w_r
          } else {
            if (method=="p") { # parametric bootstrap assuming a normal distribution
              resids <- rnorm(N,0,sd)
            } else {# non-parametric bootstrap for residuals
              std.resid = calc.StdResid(Res)$std.resid
              std.resids = sample(std.resid,replace=TRUE)
              resids = std.resids*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$input$p0 = Res$opt$par
          res.b <- try(do.call(fit.SR, res.b$input))
          if (class(res.b) != "try-error") break
        }

      return(res.b)
    })
  } else {
    # fit.SRregime
    N <- length(Res$input$SRdata$SSB)
    tmp = calc.StdResid(Res)
    sd = tmp$sigma
    std.resid = tmp$std.resid
    merged = full_join(Res$regime_pars,mutate(Res$regime_resid,Year=Res$input$SRdata$year),by="regime") %>%
      arrange(Year)
    RES = lapply(1:n, function(j){
      for (k in 1:100) {
        if (method=="d") { #data bootstrap
          Year0 = Res$input$SRdata$year[Res$input$w==1]
          Year_r = sample(Year0,replace=TRUE)
          w_r = purrr::map_dbl(Year0, function(x) sum(Year_r==x))
          res.b <- Res
          res.b$input$w <- res.b$input$SRdata$weight <- w_r
        } else {
          if (method=="p") { # parametric bootstrap assuming a normal distribution
            std.resids = rnorm(N,0,1)
          } else {# non-parametric bootstrap for residuals
            std.resids = sample(std.resid,replace=TRUE)
          }
          resids <- std.resids*sd
          pred <- obs <- resid0 <- numeric(N)
          ssb <- Res$input$SRdata$SSB
          for(i in 1:N){
            pred[i] <- SRF(ssb[i],merged$a[i],merged$b[i])
            if (i==1) {
              obs[i] <- pred[i]*exp(resids[i])
            } else {
              obs[i] <- pred[i]*exp(0*resid0[i-1])*exp(resids[i])
            }
            resid0[i] <- log(obs[i]/pred[i])
          }
          res.b <- Res
          res.b$input$SRdata$R <- obs
        }
        res.b$input$p0 = Res$opt$par
        res.b <- try(do.call(fit.SRregime, res.b$input))
        if (class(res.b) != "try-error") break
      }

      return(res.b)
    })
  }
  RES$input = arglist
  return(RES)
}

#'  profile likelihood
#' @param Res 再生産関係(\code{fit.SR()})の結果オブジェクト
#' @encoding UTF-8
#' @export
#'

prof.lik <- function(Res,a=Res$pars$a,b=Res$pars$b,sd=Res$pars$sd,rho=ifelse(Res$input$out.AR,0,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))
  validate_sr(SR = SR)
  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))
}

#' レジーム分けを考慮した再生産関係の推定
#'
#' レジームシフトが生じた年やレジームであるパラメータが共通する場合やレジームのパターンがA->B->CなのかA->B->Aなのか等が検討できる
#' @param SRdata \code{get.SRdata}で作成した再生産関係データ
#' @param SR 再生産関係 (\code{"HS"}: Hockey-stick, \code{"BH"}: Beverton-Holt, \code{"RI"}: Ricker)
#' @param method 最適化法(\code{"L2"}: 最小二乗法, \code{"L1"}: 最小絶対値法)
#' @param regime.year レジームが変わる年
#' @param regime.key レジームのパターンを表す(\code{0:2}だとA->B->Cで、\code{c(0,1,0)}だとA->B->Aのようなパターンとなる)
#' @param regime.par レジームによって変化するパラメータ(\code{c("a","b","sd")}の中から選ぶ)
#' @param length 初期値を決める際のgridの長さ
#' @param p0 \code{optim}で設定する初期値
#' @param bio_par data.frame(waa=c(100,200,300),maa=c(0,1,1),M=c(0.3,0.3,0.3)) のような生物パラメータをあらわすデータフレーム。waaは資源量を計算するときのweight at age, maaはmaturity at age, Mは自然死亡係数。これを与えると、steepnessやR0も計算して返す
#' @param plus_group hなどを計算するときに、プラスグループを考慮するか
#' @inheritParams fit.SR
#' @encoding UTF-8
#' @examples
#' \dontrun{
#' data(res_vpa)
#' SRdata <- get.SRdata(res_vpa)
#' resSRregime <- fit.SRregime(SRdata, SR="HS", method="L2",
#'                             regime.year=c(1995,2005), regime.key=c(0,1,0),
#'                             regime.par = c("a","b","sd")[2:3])
#' resSRregime$regime_pars
#' }
#' @return 以下の要素からなるリスト
#' \describe{
#' \item{\code{input}}{使用した引数のリスト}
#' \item{\code{pars}}{推定されたパラメータ}
#' \item{\code{opt}}{\code{optim}の結果オブジェクト}
#' \item{\code{resid}}{再生産関係から予測値からの加入量の残差}
#' \item{\code{loglik}}{対数尤度}
#' \item{\code{k}}{推定したパラメータ数}
#' \item{\code{AIC}}{AIC}
#' \item{\code{AICc}}{AICc}
#' \item{\code{BIC}}{BIC}
#' \item{\code{regime_pars}}{レジームごとの推定パラメータ}
#' \item{\code{regime_resid}}{レジームごとの残差}
#' \item{\code{pred}}{レジームごとの各親魚量に対する加入量の予測値}
#' \item{\code{pred_to_obs}}{観測値に対する予測値}
#' \item{\code{summary_tbl}}{観測値と予測値を合わせた表}
#' \item{\code{steepness}}{bio_parを与えたときに、steepness (h) やR0(漁獲がない場合の平均加入尾数), SB0(漁獲がない場合の平均親魚量)なども追加的に返す}
#' }
#'
#' @export
fit.SRregime <- function(
  SRdata,
  SR = "HS",
  method = "L2",
  regime.year = NULL,
  regime.key = 0:length(regime.year),
  regime.par = c("a","b","sd"),
  # sd.obs = TRUE,
  use.fit.SR = TRUE,
  length=10,  # parameter (a,b) の初期値を決めるときにgrid searchする数
  p0 = NULL,  # 初期値
  w = NULL,
  max.ssb.pred = 1.3,
  hessian = FALSE,
  bio_par = NULL,
  plus_group = TRUE,
  gamma = 0.001
) {
  argname <- ls()
  arglist <- lapply(argname,function(xx) eval(parse(text=xx)))
  names(arglist) <- argname

  #  if(!is.null(SRdata$weight)) w <- SRdata$w
  #  if(is.null(w)) w <- rep(1,length(SRdata$R))
  tmp <- check_consistent_w(w, SRdata)
  SRdata <- tmp$SRdata
  w <- arglist$w <- tmp$w

  rec <- SRdata$R
  ssb <- SRdata$SSB
  N <- length(rec)

  regime.key0 = regime.key
  unique.key = unique(regime.key)
  regime.key = sapply(1:length(regime.key0), function(i) which(unique.key == regime.key0[i]))
  regime <- a_key <- b_key <- sd_key <- rep(1,N)

  if (!is.null(regime.year)) {
    for(i in 1:length(regime.year)) {
      regime[SRdata$year>=regime.year[i]] <- regime.key[i+1]
    }
  }

  if ("a" %in% regime.par) a_key <- regime
  if ("b" %in% regime.par) b_key <- regime
  if ("sd" %in% regime.par) sd_key <- regime

  validate_sr(SR = SR, method = method)
  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)
  if (SR=="Mesnil") SRF <- function(x,a,b) 0.5*a*(x+sqrt(b^2+gamma^2/4)-sqrt((x-b)^2+gamma^2/4))

  obj.f <- function(a,b,out="nll"){ #a,bはベクトル
    resid <- NULL
    for(i in 1:N) {
      pred_rec <- SRF(ssb[i],a[a_key[i]],b[b_key[i]])
      resid <- c(resid, log(rec[i]/pred_rec))
    }
    if (method == "L2") {
      tbl = tibble(resid=resid,sd_key=sd_key,w=w) %>%
        mutate(resid2 = resid^2) %>%
        group_by(sd_key) %>%
        summarise(rss = sum(w*resid2), n = sum(w),.groups="drop") %>%
        mutate(sd = sqrt(rss/n))
    } else {
      tbl = tibble(resid=resid,sd_key=sd_key,w=w) %>%
        # mutate(resid2 = resid^2) %>%
        group_by(sd_key) %>%
        summarise(rss = sum(w*abs(resid)), n = sum(w),.groups="drop") %>%
        mutate(sd = rss/n)
    }
    SD = tbl$sd %>% as.numeric()
    nll <-
      ifelse(method=="L2",
             -sum(sapply(1:N, function(i) w[i]*dnorm(resid[i],0,SD[sd_key[i]], log=TRUE))),
             -sum(sapply(1:N, function(i) w[i]*(-log(2*SD[sd_key[i]])-abs(resid[i]/SD[sd_key[i]]))))
      )
    if (out=="nll") return(nll)
    if (out=="resid") return(resid)
    if (out=="sd") return(SD)
  }

  a_grid <- NULL
  for(i in unique(a_key)){
    a_range<-range(rec[a_key==i]/ssb[a_key==i])
    a_grid <- cbind(a_grid,seq(a_range[1],a_range[2],length=length))
  }
  b_grid <- NULL
  for(i in unique(b_key)){
    if (SR=="HS" | SR=="Mesnil") {
      b_range <- range(ssb[b_key==i])
    } else {b_range <- range(1/ssb[b_key==i])}
    b_grid <- cbind(b_grid,seq(b_range[1],b_range[2],length=length))
  }
  ab_grid <- expand.grid(data.frame(a_grid,b_grid)) %>% as.matrix()
  b_range <- apply(b_grid,2,range)

  if (is.null(p0)) {
    if (use.fit.SR) {
      fit_SR_res = fit.SR(SRdata, SR = SR, method = method, w = w, AR = 0)
      init <- c(rep(fit_SR_res$opt$par[1],max(a_key)),rep(fit_SR_res$opt$par[2],max(b_key)))
    } else {
      init_list <- sapply(1:nrow(ab_grid), function(j) {
        obj.f(a=ab_grid[j,1:max(a_key)],b=ab_grid[j,(1+max(a_key)):(max(a_key)+max(b_key))])
      })

      ab_init <- as.numeric(ab_grid[which.min(init_list),])
      init <- log(ab_init[1:max(a_key)])
      if (SR=="HS" | SR=="Mesnil") {
        for(i in unique(b_key)) {
          init <- c(init,-log(max(0.000001,(max(ssb[b_key==i])-min(ssb[b_key==i]))/max(0.000001,(ab_init[max(a_key)+i]-min(ssb[b_key==i])))-1)))
        }
      } else {
        init <- c(init, log(ab_init[(max(a_key)+1):(max(a_key)+max(b_key))]))
      }
    }
  } else {
    init <- p0
  }

  if (SR=="HS" | SR=="Mesnil") {
    obj.f2 <- function(x) {
      a <- exp(x[1:max(a_key)])
      b <- b_range[1,]+(b_range[2,]-b_range[1,])/(1+exp(-x[(max(a_key)+1):(max(a_key)+max(b_key))]))
      return(obj.f(a,b))
    }
  } else {
    obj.f2 <- function(x) obj.f(a=exp(x[1:max(a_key)]),b=exp(x[(1+max(a_key)):(max(a_key)+max(b_key))]))
  }

  opt <- optim(init,obj.f2)
  for (i in 1:100) {
    opt2 <- optim(opt$par,obj.f2)
    if (abs(opt$value-opt2$value)<1e-6) break
    opt <- opt2
  }
  opt <- optim(opt$par,obj.f2,method="BFGS",hessian=hessian)

  Res <- list()
  Res$input <- arglist
  Res$opt <- opt

  a <- exp(opt$par[1:max(a_key)])
  if (SR=="HS" | SR=="Mesnil") {
    b <- b_range[1,]+(b_range[2,]-b_range[1,])/(1+exp(-opt$par[(1+max(a_key)):(max(a_key)+max(b_key))]))
  } else {
    b <- exp(opt$par[(1+max(a_key)):(max(a_key)+max(b_key))])
  }
  sd.pred <- sd <- obj.f(a,b,out="sd")
  if (method=="L1") {
    # L1の場合sd.predとsdは定義が異なる. sdの値は825-834行目は上書きされることに注意
    sd.pred <- sd <- sqrt(2)*sd.pred
  }
  # sd <- sqrt(sum(w*obj.f(a,b,out="resid")^2)/sum(w))
  resid <- obj.f(a,b,out="resid")

  Res$obj.f <- obj.f
  Res$obj.f2 <- obj.f2
  Res$resid <- resid
  Res$pars$a <- a
  Res$pars$b <- b
  Res$pars$sd <- sd

  Res$loglik <- loglik <- -opt$value
  Res$k <- k <- length(opt$par)+max(sd_key)
  Res$AIC <- -2*loglik+2*k
  Res$AICc <- Res$AIC+2*k*(k+1)/(sum(w>0)-k-1)
  Res$BIC <- -2*loglik+k*log(sum(w>0))

  Res$regime_pars <- tibble(regime=regime.key0[regime],a=a[a_key],b=b[b_key],sd=sd[sd_key]) %>% distinct()
  Res$regime_resid <- tibble(regime=regime.key0[regime],resid = resid)
  # if (sd.obs) {
  Res$sd.pred = sd.pred
  if (method=="L1") {
      if ("sd" %in% regime.par) {
        tmp = Res$regime_resid %>% mutate(w = w, squared_resid = w*resid^2) %>%
          group_by(regime) %>% summarise(n = sum(w), RSS = sum(squared_resid)) %>%
          mutate(RMSE = sqrt(RSS/n))
        sd = as.numeric(tmp$RMSE)
      } else {
        sd = sqrt(sum(resid[w==1]^2)/sum(w))
      }
      Res$pars$sd = sd
      Res$regime_pars$sd = sd
    }
  # }
  ssb.tmp <- seq(from=0,to=max(ssb)*max.ssb.pred,length=100)
  ab_unique <- unique(cbind(a_key,b_key))
  summary_tbl = tibble(Year = SRdata$year,SSB=ssb, R = rec, Regime=regime.key0[regime], Category = "Obs")
  for (i in 1:nrow(ab_unique)) {
    R.tmp <- sapply(1:length(ssb.tmp), function(j) SRF(ssb.tmp[j],a[ab_unique[i,1]],b[ab_unique[i,2]]))
    summary_tbl = bind_rows(summary_tbl,tibble(Year=NA,SSB=ssb.tmp, R=R.tmp, Regime=unique.key[i], Category="Pred"))
  }
  summary_tbl = summary_tbl %>% mutate(Regime = factor(Regime))
  pred = dplyr::filter(summary_tbl, Category == "Pred") %>%
    dplyr::select(-Year, -Category) %>%
    dplyr::select(Regime,SSB,R)
  Res$pred <- pred
  pred_to_obs = dplyr::filter(summary_tbl, Category == "Obs") %>%
    dplyr::select(-Category) %>%
    mutate(resid = resid) %>%
    mutate(Pred = exp(log(R)-resid)) %>%
    dplyr::select(Year,SSB,R,Regime,Pred,resid)
  Res$pred_to_obs <- pred_to_obs
  Res$summary_tbl
  class(Res) <- "fit.SRregime"

  if(!is.null(bio_par)){
    if(SR!="Mesnil") par.matrix <- Res$regime_pars[c("a","b")]
    else par.matrix <- tibble(Res$regime_pars[c("a","b")],gamma=gamma)
      Res$steepness <- purrr::map_dfr(seq_len(nrow(par.matrix)),
                                function(i){
                                    calc_steepness(SR=SR,rec_pars=par.matrix[i,],M=bio_par$M,waa=bio_par$waa,maa=bio_par$maa,
                                                   plus_group=plus_group) %>%
                                        mutate(regime=Res$regime_pars$regime[i])
                                })
  }

  Res$pars <- NULL
  return(Res)
}

#' 再生産関係の推定における標準化残差を計算する関数
#' @import rmutil
#' @param resSR \code{fit.SR}か\code{fit.SRregime}のオブジェクト
#' @encoding UTF-8
#' @export
calc.StdResid = function(resSR) {
  validate_sr(res_sr = resSR)
  if(class(resSR) == "fit.SR") { #fit.SR
    if (resSR$input$method == "L2") {
      sigma = rep(sqrt(sum(resSR$resid^2)/length(resSR$resid)),length(resSR$resid))
    } else {
      sigma = rep(sqrt(2)*sum(abs(resSR$resid))/length(resSR$resid),length(resSR$resid))
    }
    sigma2 = c(sqrt(resSR$pars$sd^2/(1-resSR$pars$rho^2)), rep(resSR$pars$sd,length(resSR$resid)-1))
    std.resid = resSR$resid/sigma
    std.resid2 = resSR$resid2/sigma2
    if (resSR$input$method == "L2") {
      cumulative.prob = pnorm(std.resid,0,1)
    } else {
      cumulative.prob = rmutil::plaplace(std.resid,0,s=1/sqrt(2))
    }
    if (resSR$input$AR) {
      cumulative.prob2 = pnorm(std.resid2,0,1)
    } else {
      cumulative.prob2 = cumulative.prob
    }
    RES = tibble(sigma,sigma2,std.resid,std.resid2,cumulative.prob,cumulative.prob2)
  } else{ #fit.SRregime
    RES = dplyr::full_join(resSR$regime_pars,
                           resSR$regime_resid %>% mutate(Year = resSR$input$SRdata$year),by="regime") %>%
      dplyr::arrange(Year) %>%
      dplyr::mutate(std.resid = resid/sd) %>%
      dplyr::select(sd,std.resid) %>%
      rename(sigma=sd)
    if (resSR$input$method == "L2") {
      cumulative.prob = pnorm(RES$std.resid,0,1)
    } else {
      cumulative.prob = rmutil::plaplace(RES$std.resid,0,s=1/sqrt(2))
    }
    RES = RES %>% mutate(cumulative.prob=cumulative.prob)
  }
  return(RES)
}


#' 再生産関係の残差の確率分布に関するチェック図を出力する関数
#'
#' 1) 正規性をチェックするための検定結果・ヒストグラム、2) 累積確率分布のヒストグラム、3) 一様分布を使用した累積確率分布のQQプロttの3つが出力される
#' 標準化した残差を使用する
#' 累積確率分布はL2の場合は正規分布、L1の場合はラプラス分布を使用する
#' 自己相関の外側推定の場合は、2段階で推定しているため、ファイルが2つ出力される
#' @import EnvStats
#' @import rmutil
#' @inheritParams calc.StdResid
#' @param resSR \code{fit.SR}か\code{fit.SRregime}のオブジェクト
#' @param output pngファイルに出力するか否か
#' @encoding UTF-8
#' @export
check.SRdist = function(resSR,test.ks=TRUE,output=FALSE,filename = "SR_error_dist") {
  std.resid.res = calc.StdResid(resSR)
  main_name=paste0(resSR$input$SR," ",resSR$input$method," ")

  if (class(resSR)=="fit.SR" && resSR$input$AR && isTRUE(resSR$input$out.AR)) {
    validate_sr(res_sr = resSR)
    for(i in 1:2) {
      if (i==1) {
        std.resid = std.resid.res$std.resid
        cumulative.prob = std.resid.res$cumulative.prob
        main_name2 = "Std. Deviance to SR"
      } else {
        std.resid = std.resid.res$std.resid2
        cumulative.prob = std.resid.res$cumulative.prob2
        main_name2 = "Standard. Resid."
      }
      # pdf(file = paste0(filename,"(",main_name2,").pdf"), width=15,height=5)
      if (output) png(file = paste0(filename,"(",main_name2,").png"), width=15, height=5, res=432, units='in')
      par(pch=1,lwd = 2, mfrow=c(1,3),cex=1)
      check1 <- shapiro.test(std.resid)
      hist(std.resid,main = paste0(main_name,main_name2),
           xlab = "Standardized Deviance",freq=FALSE,col="gray")
      X <- seq(min(std.resid)*1.3,max(std.resid)*1.3,length=200)
      points(X,dnorm(X,0,1),col=2,lwd=3,type="l")
      if (resSR$input$method=="L1" && i==1) {
        points(X,rmutil::dlaplace(X,0,1/sqrt(2)),col="blue",lwd=2,type="l",lty=2)
      }
      mtext(text=" P value",adj=1,line=-1,lwd=2,font=2)
      mtext(text=sprintf(" SW: %1.3f",check1$p.value),adj=1,line=-2)
      if (test.ks) {
        if (i==1 && resSR$input$method=="L1") {
          check2 <- ks.test(std.resid,y=function(x) rmutil::plaplace(x,m=0,s=1/sqrt(2)))
        } else {
          check2 <- ks.test(std.resid,y="pnorm",sd=1)
        }
        mtext(text=sprintf(" KS: %1.3f",check2$p.value),adj=1,line=-3)
      }
      hist(cumulative.prob, xlab="Cumulative probability",
           main = paste0(main_name,"Cumlative Prob. Dist."),freq=FALSE,breaks=seq(0,1,by=0.1),col="gray")
      # if (test.ks) {
      #   check2 <- ks.test(cumulative.prob,y="punif",min=0,max=1)
      #   mtext(text=" P value",adj=1,line=-1,lwd=2,font=2)
      #   mtext(text=sprintf(" KS: %1.3f",check2$p.value),adj=1,line=-2)
      # }
      EnvStats::qqPlot(x=cumulative.prob,distribution="unif", param.list = list(min = 0, max = 1),
                       add.line = TRUE,qq.line.type = "0-1",line.lwd=1.5,cex=1.2,line.col="red",
                       main = paste0(main_name,"Uniform QQ plot"),ylab="Quatiles of Cumlative.Prob.",
                       xlab = "Quantiles of Uniform Dist.")
      if (output) dev.off()
    }
  } else {
    if (is.null(std.resid.res$std.resid2)) {
      std.resid = std.resid.res$std.resid
      cumulative.prob = std.resid.res$cumulative.prob
    } else {
      std.resid = std.resid.res$std.resid2
      cumulative.prob = std.resid.res$cumulative.prob2
    }
    # pdf(file = paste0(filename,".pdf"), width=15,height=5)
    if (output) png(file = paste0(filename,".png"), width=15, height=5, res=432, units='in')
    par(pch=1,lwd = 2, mfrow=c(1,3),cex=1)
    check1 <- shapiro.test(std.resid)
    hist(std.resid,main = paste0(main_name,"Standard. Resid."),
         xlab = "Standardized Residual",freq=FALSE,col="gray")
    X <- seq(min(std.resid)*1.3,max(std.resid)*1.3,length=200)
    points(X,dnorm(X,0,1),col=2,lwd=3,type="l")
    if (resSR$input$method=="L1") {
      points(X,rmutil::dlaplace(X,0,1/sqrt(2)),col="blue",lwd=2,type="l",lty=2)
    }
    mtext(text=" P value",adj=1,line=-1,lwd=2,font=2)
    mtext(text=sprintf(" SW: %1.3f",check1$p.value),adj=1,line=-2)
    if (test.ks) {
      if (resSR$input$method=="L1"){
        check2 <- ks.test(std.resid,y=function(x) rmutil::plaplace(x,m=0,s=1/sqrt(2)))
      } else {
        check2 <- ks.test(std.resid,y="pnorm",mean=0,sd=1)
      }
      mtext(text=sprintf(" KS: %1.3f",check2$p.value),adj=1,line=-3)
    }
    hist(cumulative.prob, xlab="Cumulative probability",
         main = paste0(main_name,"Cumlative Prob. Dist."),freq=FALSE,breaks=seq(0,1,by=0.1),col="gray")
    # if (test.ks) {
    #   check2 <- ks.test(cumulative.prob,y="punif",min=0,max=1)
    #   mtext(text=" P value",adj=1,line=-1,lwd=2,font=2)
    #   mtext(text=sprintf(" KS: %1.3f",check2$p.value),adj=1,line=-2)
    # }
    EnvStats::qqPlot(x=cumulative.prob,distribution="unif", param.list = list(min = 0, max = 1),
                     add.line = TRUE,qq.line.type = "0-1",line.lwd=1.5,cex=1.2,line.col="red",
                     main = paste0(main_name,"Uniform QQ plot"),ylab="Quatiles of Cumlative.Prob.",
                     xlab = "Quantiles of Uniform Dist.")
    if (output) dev.off()
  }
}

#' 再生産関係の残差から事後的に1次の自己相関係数を推定する関数
#'
#' 計算は\code{fit.SR}の\code{out.AR}と同じであるが、AICcも見れるほか、\code{fit.SRregime}にも対応している
#' \code{fit.SRregime}の場合、各レジームの最初の年は初期値となる(つまり前のレジームの最後の年からの残差を引きずらない)
#' @param resSR \code{fit.SR}か\code{fit.SRregime}のオブジェクト
#' @param per_regime 自己相関係数をレジームごとに推定するか (\code{TRUE}) 否か
#' @encoding UTF-8
#' @examples
#' \dontrun{
#' data(res_vpa)
#' SRdata <- get.SRdata(res_vpa)
#' data(res_vpa)
#' SRdata <- get.SRdata(res_vpa)
#' resSR <- fit.SR(SRdata, SR = c("HS","BH","RI")[1],
#'                 method = c("L1","L2")[1], AR = 0,
#'                 out.AR = FALSE)
#' resSR_post = calc.residAR(resSR)
#' resSR_post$AICc
#' resSR_post$pars
#' resSRregime <- fit.SRregime(SRdata, SR="HS", method="L2",
#'                             regime.year=c(1995,2005), regime.key=c(0,1,0),
#'                             regime.par = c("a","b","sd")[2:3])
#' resSRregime_post = calc.residAR(resSRregime, per_regime = TRUE)
#' resSRregime_post$AICc
#' resSRregime_post$regime$pars
#' }
#'
#' @export
calc.residAR = function(resSR, per_regime=TRUE, output=TRUE, filename="residARouter") {
  RES = list()
  if (class(resSR) == "fit.SR") { #fit.SR
    validate_sr(res_sr = resSR)
    if (resSR$input$AR && !isTRUE(resSR$input$out.AR)) {
      warning("This function is meaningless when AR=TRUE & out.AR=FALSE")
    }
    deviance = resSR$resid
    arimares0 = arima(deviance,method="ML",include.mean=FALSE,order=c(0,0,0))
    arimares1 = arima(deviance,method="ML",include.mean=FALSE,order=c(1,0,0))
    nn = length(deviance)
    k = c("AR(0)"=1,"AR(1)"=2)
    loglik = c("AR(0)"=arimares0$loglik,"AR(1)"=arimares1$loglik)
    AIC = c("AR(0)"=arimares0$aic,"AR(1)"=arimares1$aic)
    AICc = AIC + 2*k*(k+1)/(nn-k-1)
    BIC = AIC + k*log(nn)
    rho = as.numeric(arimares1$coef)
    sd = sqrt(arimares1$sigma2)
    resid2 = arimares1$residuals
    RES$pars = resSR$pars
    RES$pars$rho = rho
    RES$pars$sd = sd
    RES$arima0 = arimares0
    RES$arima1 = arimares1
    RES$resid = resSR$resid
    RES$resid2 = resid2
    RES$loglik = loglik
    RES$k = k
    RES$AIC = AIC
    RES$AICc = AICc
    RES$BIC = BIC
  } else {  #fit.SRregime
    tbl = resSR$regime_pars %>% full_join(resSR$regime_resid %>% mutate(Year = resSR$input$SRdata$year)) %>%
      arrange(Year) %>% mutate(Shift = if_else(Year %in% c(min(Year), resSR$input$regime.year),1,0))
    tbl = tbl %>%
      mutate(regime_id = purrr::map_dbl(regime, function(x) which(x==unique(resSR$regime_pars$regime))))
    deviance = tbl$resid; shift = tbl$Shift; regime_id = tbl$regime_id;
    w = resSR$input$w;nn = nrow(tbl);
    obj.f = function(x,out="nll") {
      rho = (exp(x)-1.0)/(exp(x)+1.0);resid2 = c();
      if (per_regime) {
        for (i in 1:nn) {
          if (shift[i]==1) {
            resid2 = c(resid2,sqrt(1-rho[regime_id[i]]^2)*deviance[i])
          } else {
            resid2 = c(resid2,deviance[i]-rho[regime_id[i-1]]*deviance[i-1])
          }
        }
      } else {
        for (i in 1:nn) {
          if (shift[i]==1) {
            resid2 = c(resid2,sqrt(1-rho^2)*deviance[i])
          } else {
            resid2 = c(resid2,deviance[i]-rho*deviance[i-1])
          }
        }
      }
      RSS = tbl %>% mutate(resid2 = resid2,w=w) %>% group_by(sd) %>%
        summarise(var=sum(w*resid2^2),n=sum(w)) %>%
        mutate(sigma = sqrt(var/n))
      tbl2 = left_join(tbl,RSS,by="sd")
      sigma = tbl2$sigma
      if (out=="nll") {
        return(-sum(dnorm(0,resid2,sigma,log=TRUE))) } else {
          if (out=="resid") {
            return(resid2) } else {
              return(sigma)
            }
        }
    }
    if (per_regime) {
      opt = optim(rep(0,max(tbl$regime_id)),obj.f)
      for (i in 1:100) {
        opt2 <- optim(opt$par,obj.f)
        if (abs(opt$value-opt2$value)<1e-6) break
        opt <- opt2
      }
      opt <- optim(opt$par,obj.f,method="BFGS")
    } else {
      opt = optim(0,obj.f,method="Brent",lower=-20,upper=20)
    }
    rho = (exp(opt$par)-1)/(exp(opt$par)+1)
    loglik1 = -opt$value
    resid2 = obj.f(x=opt$par,out="resid")
    sigma = obj.f(x=opt$par,out="sigma")
    tbl3 = tbl %>% mutate(resid2 = resid2,sigma=sigma)
    RES$pars = resSR$pars
    RES$pars$rho = rho
    if (per_regime) rho = rho[tbl3$regime_id]
    tbl3 = tbl3 %>% mutate(rho = rho)
    RES$regime_pars = tbl3 %>% dplyr::select(regime,a,b,sigma,rho) %>% distinct() %>% rename(sd = sigma)
    RES$regime_resid = tbl3 %>% dplyr::select(regime,resid,resid2)
    RES$pars$sd = unique(RES$regime_pars$sd)
    RES$opt = opt
    if (per_regime) {
      loglik0 = -obj.f(rep(0,max(tbl$regime_id)))
    } else {
      loglik0 = -obj.f(0)
    }
    k = c("AR(0)"=length(resSR$pars$sd),"AR(1)"=length(resSR$pars$sd)+length(opt$par))
    loglik = c("AR(0)"=loglik0,"AR(1)"=loglik1)
    AIC = -2*loglik+2*k
    AICc = AIC + 2*k*(k+1)/(nn-k-1)
    BIC = AIC + k*log(nn)
    RES$loglik = loglik
    RES$k = k
    RES$AIC = AIC
    RES$AICc = AICc
    RES$BIC = BIC
  }
  if (output) {
    capture.output(RES, file = paste0(filename,".txt"))
  }
  return(RES)
}

#' 再生産関係における残差の時系列の自己相関等についてのプロット
#'
#' 1) 残差のトレンド、2) \code{acf}関数による自己相関係数のプロット、3) Ljung-Box検定におけるP値の3つの図を出力
#' @inheritParams calc.StdResid
#' @param resSR \code{fit.SR}か\code{fit.SRregime}のオブジェクト
#' @param use.resid 再生産関係との残差 (deviance \code{resid}) を使うか (1:default)、自己相関を除いた残差 (\code{resid2}) を使うか (2)
#' @param output pngファイルに出力するか否か
#' @encoding UTF-8
#' @export
autocor.plot = function(resSR,use.resid=1,lag.max=NULL,output = FALSE,filename = "Residual_trend",pch=16,lwd=2,cex=1.2,cex.main=1.2,cex.lab=1.2,...){
  Year = resSR$input$SRdata$year
  if (output) png(file = paste0(filename,".png"), width=15, height=5, res=432, units='in')
  par(pch=pch,lwd = lwd, mfrow=c(1,3),cex=cex)
  if (class(resSR) == "fit.SR") { #fit.SR
    validate_sr(res_sr = resSR)
    if (use.resid==1) {
      Resid = resSR$resid
      plot(Year,Resid,pch=pch,main="",xlab="Year",ylab="Deviance",cex.lab=cex.lab,...)
      title("Time series of deviance to SR",cex.main=cex.main)
    } else {
      Resid = resSR$resid2
      Resid[1] = sqrt(1-resSR$pars$rho^2)*Resid[1]
      plot(Year,Resid,pch=pch,main="",xlab="Year",ylab="Residual",cex.lab=cex.lab,...)
      title("Time series of Residuals",cex.main=cex.main)
    }
  } else { #fit.SRregime
    message("Standardized residuals are used for 'fit.SRregime'")
    table = calc.StdResid(resSR)
    Resid = table$std.resid
    plot(Year,Resid,pch=pch,main="",xlab="Year",ylab="Std.Residual",cex.lab=cex.lab,...)
    title("Time series of Standardized Residuals",cex.main=cex.main)
    abline(v=c(resSR$input$regime.year)-0.5,lty=3,col="blue")
  }
  abline(0,0,lty=2)
  par(new=T)
  scatter.smooth(Year, Resid, lpars=list(col="red",lwd=lwd),ann=F,axes=FALSE)

  if (is.null(lag.max)) lag.max = 10*log10(length(Resid))
  ac.res <- acf(Resid,plot=FALSE,lag.max=lag.max)
  plot(ac.res,main="",lwd=lwd,cex=cex,cex.lab=cex.lab,...)
  title("Autocorrelation (rho vs. lag)",cex.main=cex.main)

  p = c()
  col = c()
  for(i in 1:lag.max){
    LBtest = Box.test(Resid,lag = i,type="Ljung")
    p = c(p,LBtest$p.value)
    col = c(col,isTRUE(LBtest$p.value<0.05)+1)
  }
  plot(p,pch=pch,...,ylim=c(0,max(c(p,0.2))),col=col,xlab="Lag",ylab="P value",cex.lab=cex.lab,...)
  abline(0.05,0.,lty=2,col="blue")
  title("Ljung-Box test",cex.main=cex.main)
  if (output) dev.off()
}

#' 再生産関係の残差ブートストラップをプロットする関数
#'
#' ブートストラップによる各パラメータのヒストグラムと、それらパラメータをつかった再生産関係が描画される
#'
#' ヒストグラムの描画ではとその中央値、上側・下側CIが(デフォルトのggplt=Tでは加えて平均値が)推定値とともに表示される
#'
#' ヒストグラムのプロットにはggplotをつかうオプションggplt=TRUEをデフォルトにしているが、ggplotを使わない作図もggplt=Fで実行可能
#' 
#' ggplotで描画する場合でRstudio利用時にはplotをZoomで表示にしないと描画されないことがあるので注意
#' @param boot.res \code{boot.SR}のオブジェクト
#' @param CI プロットする信頼区間
#' @param output pngファイルに出力するか否か
#' @param filename ファイル名
#' @encoding UTF-8
#' @export
bootSR.plot = function(boot.res, CI = 0.8,output = FALSE,filename = "boot",lwd=1.2,pch=1,ggplt=TRUE,...) {
  res_base = boot.res$input$Res
  if (class(boot.res$input$Res)=="fit.SR") {
    validate_sr(res_sr = boot.res$input$Res)
    # for fit.SR ----
    # parameter histogram ----
    if(ggplt){ # ggplot (plot histograms of a,b,B0,h,R0,rho,SB0,sd)
      if(is.null(convert_SR_tibble(boot.res[[1]]))) print("Do fit.SR with argument(bio.par) to calculate steepness.")
      N <- boot.res$input$n
      res_boot_par_tibble <- purrr::map_dfr(boot.res[1:N], convert_SR_tibble) %>% dplyr::filter(type=="parameter"&name!="SPR0")
      # AR=0のときrhoのプロットを除く
      if(boot.res$input$Res$input$AR==0) res_boot_par_tibble <- res_boot_par_tibble %>% dplyr::filter(name!="rho")

      res_boot_par_tibble_summary <- res_boot_par_tibble %>%  group_by(name) %>% summarise(median=median(value),mean=mean(value),CI10=quantile(value,0.1),CI90=quantile(value,0.9)) %>% mutate(name_with_CI = str_c(name," (",round(CI10,2),"-",round(CI90,2),")")) %>% pivot_longer(cols=-c(name,name_with_CI), names_to="stats")

      res_boot_par_base <- res_boot_par_tibble_summary %>% filter(stats=="median") #推定結果のために中央値から形だけ持ってくる
      res_fitSRtibble <- convert_SR_tibble(boot.res[["input"]]$Res) %>% dplyr::filter(type=="parameter"&name!="SPR0")


      for(i in 1:length(res_boot_par_base)){
        for(j in 1:length(res_boot_par_base)){
          if(res_boot_par_base$name[i] ==res_fitSRtibble$name[j])
            res_boot_par_base$value[i] <- res_fitSRtibble$value[j]
        }
      }
      res_boot_par_base$stats <- "estimated"

      res_boot_par_tibble <- res_boot_par_tibble_summary %>% select(name, name_with_CI) %>% right_join(res_boot_par_tibble)

      plot_col <- c("blue","blue","red","darkgreen","green")
      base_linetype <- c("solid")
      bootsr_linetype <- c("dashed","dashed","longdashed","solid")
      plot_bootsr_linetype <- rep(bootsr_linetype,length(levels(as.factor(res_boot_par_tibble$name))))
      #legend.labels <- c("estimated", "CI10", "CI90","mean","median")

      if (boot.res$input$method=="d") {
        plot_bootsr_title<- paste0("Data Bootstrap")
      } else {
        if (boot.res$input$method=="p") {
          plot_bootsr_title<- paste0("Parametric Bootstrap")
        } else {
          plot_bootsr_title<- paste0("Non-Parametric Bootstrap")
        }
      }

      boot_par_hist <-ggplot(res_boot_par_tibble) + geom_histogram(aes(x=value)) + facet_wrap(.~fct_inorder(name_with_CI), scale="free")+theme_SH(legend.position="bottom")  + labs(title=plot_bootsr_title)+
        geom_vline(data=res_boot_par_base, mapping = aes(xintercept=value,color=stats),linetype=base_linetype) +
        geom_vline(data=res_boot_par_tibble_summary, mapping = aes(xintercept=value,color=stats),linetype="dashed") +
        scale_color_manual(name="stats",values = plot_col) #+
      #scale_linetype_manual(name="",values = plot_bootsr_linetype) #+ #scale_color_discrete(name="stats",breaks=legend.labels)

      print(boot_par_hist)

      if (output) ggsave(file = paste0(filename,"_pars.png"), plot=boot_par_hist, width=10, height=10,  units='in')

    }
    else { #plot not using ggplot
      if (output) png(file = paste0(filename,"_pars.png"), width=10, height=10, res=432, units='in')
      par(pch=pch,lwd = lwd, mfrow=c(2,2))
      jmax = ifelse(boot.res$input$Res$pars$rho==0,3,4)
      for (j in 1:jmax) {
        par0 = c("a","b","sd","rho")[j]

        hist(sapply(1:boot.res$input$n, function(i) boot.res[[i]]$pars[,par0]),xlab=par0,ylab="Frequency",main="",col="gray")
        abline(v=boot.res$input$Res$pars[,par0],col=2,lwd=3)
        abline(v=median(sapply(1:boot.res$input$n, function(i) boot.res[[i]]$pars[,par0])),col=3,lwd=3,lty=2)
        arrows(quantile(sapply(1:boot.res$input$n, function(i) boot.res[[i]]$pars[,par0]),0.5*(1-CI)),0,
               quantile(sapply(1:boot.res$input$n, function(i) boot.res[[i]]$pars[,par0]),0.5*(1-CI)+CI),0,
               col=4,lwd=3,code=3)
        legend("topright",
               legend=c("Estimate","Median","CI(0.8)"),lty=1:2,col=2:4,lwd=2,ncol=1,cex=1)
        if (boot.res$input$method=="d") {
          title(paste0(par0," in Data Bootstrap"))
        } else {
          if (boot.res$input$method=="p") {
            title(paste0(par0," in Parametric Bootstrap"))
          } else {
            title(paste0(par0," in Non-Parametric Bootstrap"))
          }
        }
      }
      if (output) dev.off()

      # steepness histogram (if available) ----
      if(!is.null(boot.res[[1]]$steepness)){
        if (output) png(file = paste0(filename,"_pars_steepness.png"), width=10, height=10, res=432, units='in')
        par(pch=pch,lwd = lwd, mfrow=c(2,2))
        for (j in 1:4) {
          par0 = c("SB0","R0","B0","h")[j]

          hist(sapply(1:boot.res$input$n, function(i) boot.res[[i]]$steepness[,par0]),xlab=par0,ylab="Frequency",main="",col="gray")
          abline(v=boot.res$input$Res$steepness[,par0],col=2,lwd=3)
          abline(v=median(sapply(1:boot.res$input$n, function(i) boot.res[[i]]$steepness[,par0])),col=3,lwd=3,lty=2)
          arrows(quantile(sapply(1:boot.res$input$n, function(i) boot.res[[i]]$steepness[,par0]),0.5*(1-CI)),0,
                 quantile(sapply(1:boot.res$input$n, function(i) boot.res[[i]]$steepness[,par0]),0.5*(1-CI)+CI),0,
                 col=4,lwd=3,code=3)
          legend("topright",
                 legend=c("Estimate","Median","CI(0.8)"),lty=1:2,col=2:4,lwd=2,ncol=1,cex=1)
          if (boot.res$input$method=="d") {
            title(paste0(par0," in Data Bootstrap"))
          } else {
            if (boot.res$input$method=="p") {
              title(paste0(par0," in Parametric Bootstrap"))
            } else {
              title(paste0(par0," in Non-Parametric Bootstrap"))
            }
          }
        }
        if (output) dev.off()
      }

    }


    # SR curve ----
    par(mfrow=c(1,1),pch=pch,lwd = lwd)
    if (output) png(file = paste0(filename,"_SRcurve.png"), width=10, height=7.5, res=432, units='in')
    data_SR = boot.res$input$Res$input$SRdata
    plot(data_SR$R ~ data_SR$SSB, cex=2, type = "p",xlab="SSB",ylab="R",pch=1,
         main="",ylim=c(0,max(data_SR$R)*1.3),xlim=c(0,max(data_SR$SSB)*1.3))
    if (boot.res$input$method=="d") {
      title("Data Bootstrap")
    } else {
      if (boot.res$input$method=="p") {
        title("Parametric Bootstrap for Residuals")
      } else {
        title("Non-Parametric Bootstrap for Residuals")
      }
    }
    points(rev(data_SR$SSB)[1],rev(data_SR$R)[1],col=1,type="p",lwd=3,pch=16,cex=2)
    for (i in 1:boot.res$input$n) {
      points(boot.res[[i]]$pred$SSB,boot.res[[i]]$pred$R,type="l",lwd=2,col=rgb(0,0,1,alpha=0.1))
    }
    points(res_base$pred$SSB,res_base$pred$R,col=2,type="l",lwd=3)
    if (output) dev.off()
  } else {
    # fit.SRregime ----

    regime_unique = boot.res$input$Res$regime_pars$regime
    obs_data = boot.res$input$Res$pred_to_obs

    if(ggplt){ # plot using ggplot
      res_base = boot.res$input$Res
      N <- boot.res$input$n
      res_boot_par_tibble <- purrr::map_dfr(boot.res[1:N], convert_SR_tibble) %>% dplyr::filter(type=="parameter"&name!="SPR0")

      legend.labels <- c("estimated", "CI10", "CI90","mean","median")
      plot_col <- c("blue","blue","red","darkgreen","green")
      base_linetype <- c("solid")
      bootsr_linetype <- c("dashed","dashed","longdashed","solid")
      plot_bootsr_linetype <- rep(bootsr_linetype,length(levels(as.factor(res_boot_par_tibble$name))))

      regime.num <- length(levels(as.factor(res_boot_par_tibble$regime)))-1

      for(k in 0:regime.num){
        if (boot.res$input$method=="d") {
          plot_bootsr_title<- paste0("Data Bootstrap regime",k)
        } else {
          if (boot.res$input$method=="p") {
            plot_bootsr_title<- paste0("Parametric Bootstrap regime",k)
          } else {
            plot_bootsr_title<- paste0("Non-Parametric Bootstrap regime",k)
          }
        }


        res_boot_par_tibble_regime<- res_boot_par_tibble %>% filter(regime==k)

        res_boot_par_tibble_summary <- res_boot_par_tibble_regime %>%  group_by(name) %>% summarise(median=median(value),mean=mean(value),CI10=quantile(value,0.1),CI90=quantile(value,0.9)) %>% mutate(name_with_CI = str_c(name," (",round(CI10,2),"-",round(CI90,2),")")) %>% pivot_longer(cols=-c(name,name_with_CI), names_to="stats")

        res_boot_par_base <- res_boot_par_tibble_summary %>% filter(stats=="median")
        res_fitSRtibble <- convert_SR_tibble(res_base) %>% dplyr::filter(type=="parameter"&name!="SPR0") %>% filter(regime==k)
        for(i in 1:nrow(res_boot_par_base)){
          for(j in 1:nrow(res_boot_par_base)){
            if(res_boot_par_base$name[i] == res_fitSRtibble$name[j])
              res_boot_par_base$value[i] <- res_fitSRtibble$value[j]
          }
        }
        res_boot_par_base$stats <- "estimated"

        res_boot_par_tibble_regime <- res_boot_par_tibble_summary %>% select(name, name_with_CI) %>% right_join(res_boot_par_tibble_regime)

        boot_par_hist <- ggplot(res_boot_par_tibble_regime) + geom_histogram(aes(x=value)) + facet_wrap(.~fct_inorder(name_with_CI), scale="free")+theme_SH(legend.position="bottom")  + labs(title=plot_bootsr_title)+
          geom_vline(data=res_boot_par_base, mapping = aes(xintercept=value,color=stats),linetype=base_linetype) +
          geom_vline(data=res_boot_par_tibble_summary, mapping = aes(xintercept=value,color=stats),linetype="dashed") +
          scale_color_manual(name="stats",values = plot_col)

        print(boot_par_hist)

        if (output) ggsave(file = paste0(filename,"_regime",k,"_pars.png"), plot=boot_par_hist, width=10, height=10,  units='in')

      }

    }
    else{ # plot not using ggplot

    # histogram ----
    if (output) png(file = paste0(filename,"_pars.png"), width=15, height=5*nrow(boot.res$input$Res$regime_pars), res=432, units='in')
    par(lwd = lwd, mfrow=c(nrow(boot.res$input$Res$regime_pars),3))
    for (ii in 1:nrow(boot.res$input$Res$regime_pars)) {
      regime = boot.res$input$Res$regime_pars$regime[ii]
      jmax = 3
      for (j in 1:jmax) {
        par0 = c("a","b","sd","rho")[j]
        boot_pars = sapply(1:boot.res$input$n, function(i) as.numeric(boot.res[[i]]$regime_pars[ii,par0]))
        hist(boot_pars,xlab=par0,ylab="Frequency",main="",col="gray")
        abline(v=as.numeric(boot.res$input$Res$regime_pars[ii,par0]),col=2,lwd=3)
        abline(v=median(boot_pars),col=3,lwd=3,lty=2)
        arrows(quantile(boot_pars,0.5*(1-CI)),0,
               quantile(boot_pars,0.5*(1-CI)+CI),0,
               col=4,lwd=3,code=3)
        legend("topright",
               legend=c("Estimate","Median","CI(0.8)"),lty=1:2,col=2:4,lwd=2,ncol=1,cex=1)
        if (boot.res$input$method=="d") {
          title(paste0(par0," of Regime",regime," in Data Bootstrap"))
        } else {
          if (boot.res$input$method=="p") {
            title(paste0(par0," of Regime",regime," in Para Bootstrap"))
          } else {
            title(paste0(par0," of Regime",regime," in Non-Para Bootstrap"))
          }
        }

      }
    }
    if (output) dev.off()

    # histogram (steepness) ----
    if(!is.null(boot.res[[1]]$steepness)){
      if (output) png(file = paste0(filename,"_pars_steepness.png"), width=15, height=5*nrow(boot.res$input$Res$regime_pars), res=432, units='in')
      par(lwd = lwd, mfrow=c(nrow(boot.res$input$Res$regime_pars),4))
      for (ii in 1:nrow(boot.res$input$Res$regime_pars)) {
        regime = boot.res$input$Res$regime_pars$regime[ii]
        jmax = 4
        for (j in 1:jmax) {
          par0 = c("SB0","R0","B0","h")[j]
          boot_pars = sapply(1:boot.res$input$n, function(i) as.numeric(boot.res[[i]]$steepness[ii,par0]))
          hist(boot_pars,xlab=par0,ylab="Frequency",main="",col="gray")
          abline(v=as.numeric(boot.res$input$Res$steepness[ii,par0]),col=2,lwd=3)
          abline(v=median(boot_pars),col=3,lwd=3,lty=2)
          arrows(quantile(boot_pars,0.5*(1-CI)),0,
                 quantile(boot_pars,0.5*(1-CI)+CI),0,
                 col=4,lwd=3,code=3)
          legend("topright",
                 legend=c("Estimate","Median","CI(0.8)"),lty=1:2,col=2:4,lwd=2,ncol=1,cex=1)
          if (boot.res$input$method=="d") {
            title(paste0(par0," of Regime",regime," in Data Bootstrap"))
          } else {
            if (boot.res$input$method=="p") {
              title(paste0(par0," of Regime",regime," in Para Bootstrap"))
            } else {
              title(paste0(par0," of Regime",regime," in Non-Para Bootstrap"))
            }
          }

        }
      }
      if (output) dev.off()
    }
    }

    # SR curve -----
    if (output) png(file = paste0(filename,"_SRcurve.png"), width=7.5*length(regime_unique), height=7.5, res=432, units='in')
    par(mfrow=c(1,length(regime_unique)))
    for(i in 1:length(regime_unique)) {
      use_data = dplyr::filter(obs_data,Regime == regime_unique[i])
      plot(use_data$R ~ use_data$SSB, cex=2, type = "p",pch=1,xlab="SSB",ylab="R",
           main="",ylim=c(0,max(use_data$R)*1.3),xlim=c(0,max(obs_data$SSB)*1.3))
      if (boot.res$input$method=="d") {
        title(paste0("Data Bootstrap for Regime ",regime_unique[i]))
      } else {
        if (boot.res$input$method=="p") {
          title(paste0("Para Bootstrap for Regime ",regime_unique[i]))
        } else {
          title(paste0("Non-para Bootstrap for Regime ",regime_unique[i]))
        }
      }
      for (j in 1:boot.res$input$n) {
        pred_data = boot.res[[j]]$pred %>% filter(Regime == regime_unique[i])
        points(pred_data$SSB,pred_data$R,type="l",lwd=2,col=rgb(0,0,1,alpha=0.1))
      }
      pred_data = res_base$pred %>% filter(Regime == regime_unique[i])
      points(pred_data$SSB,pred_data$R,col=2,type="l",lwd=3)
    }
    if (output) dev.off()
  }

}

#' 再生産関係のジャックナイフ解析
#'
#' 結果のプロットもこの関数で行う
#' @param resSR \code{fit.SR}か\code{fit.SRregime}のオブジェクト
#' @param use.p0 初期値に\code{resSR}の結果を使うか否か
#' @param is.plot プロットするかどうか
#' @param output pngファイルに出力するか否か
#' @param filename ファイル名
#' @encoding UTF-8
#' @examples
#' \dontrun{
#' data(res_vpa)
#' SRdata <- get.SRdata(res_vpa)
#' resSR <- fit.SR(SRdata, SR = c("HS","BH","RI")[1],
#'                 method = c("L1","L2")[2], AR = 1,
#'                 out.AR = TRUE)
#' res_jackSR <- jackknife.SR(resSR,output = T)
#'
#' # if calculate steepness
#' bio_par <- derive_biopar(res_obj=res_vpa,derive_year = 2010)
#' resSR <- fit.SR(SRdata, SR = c("HS","BH","RI")[1],
#'                 method = c("L1","L2")[2], AR = 1,
#'                 out.AR = TRUE,bio_par=bio_par)
#' res_jackSR <- jackknife.SR(resSR,output = T)
#' }
#' @export
#'

jackknife.SR = function(resSR,is.plot=TRUE,use.p0 = TRUE, output=FALSE,filename = "jackknife",ylim.range = c(0,1.5),pch=19,cex=1.1,...) {

  used_data <- which(resSR$input$w==1)
  if(is.null(resSR$input$SRdata$weight)) resSR$input$SRdata$weight <- resSR$input$w
  RES = lapply(used_data, function(i){
    jack <- resSR
    jack$input$w[i] <- jack$input$SRdata$weight[i] <- 0
    if (use.p0) jack$input$p0 <- resSR$opt$par
    if (class(resSR)=="fit.SR") {
      validate_sr(res_sr = resSR)
      do.call(fit.SR,jack$input)
    } else {
      if (class(resSR)=="fit.SRregime") {
        do.call(fit.SRregime,jack$input)
      } else {
        do.call(fit.SRalpha,jack$input)
      }
    }
  })
  if (is.plot) {
    jack.res <- RES
    data_SR = resSR$input$SRdata
    years <- data_SR$year[used_data]
    # no regime ----
    if (class(resSR)=="fit.SR" || class(resSR)=="fit.SRalpha") {         # plot parameters ----
      if (output) png(file = paste0(filename,"_pars.png"), width=10, height=10, res=432, units='in')
      par(mfrow=c(2,2),mar=c(3,3,2,2),oma=c(3,3,2,2),pch=pch,cex=cex)
      plot(years,sapply(1:length(used_data), function(i) jack.res[[i]]$pars$a),type="b",
           xlab="Year removed",ylab="",main="a in jackknife",ylim=resSR$pars$a*ylim.range)
      abline(resSR$pars$a,0,lwd=2,col=2,lty=2)

      plot(years,sapply(1:length(used_data), function(i) jack.res[[i]]$pars$b),type="b",
           xlab="Year removed",ylab="",main="b in jackknife",ylim=resSR$pars$b*ylim.range)
      abline(resSR$pars$b,0,lwd=2,col=2,lty=2)

      plot(years,sapply(1:length(used_data), function(i) jack.res[[i]]$pars$sd),type="b",
           xlab="Year removed",ylab="",main="sd in jackknife",ylim=resSR$pars$sd*ylim.range)
      abline(resSR$pars$sd,0,lwd=2,col=2,lty=2)

      if (class(resSR)=="fit.SR") {
        if (resSR$input$AR==1) {
          plot(years,sapply(1:length(used_data), function(i) jack.res[[i]]$pars$rho),type="b",
               xlab="Year removed",ylab="",main="rho in jackknife",ylim=resSR$pars$rho*ylim.range)
          abline(resSR$pars$rho,0,lwd=2,col=2,lty=2)
        }
      }
      if (output) dev.off()
      if(!is.null(resSR$steepness)){ #steepness----
        if (output) png(file = paste0(filename,"_steepness.png"), width=10, height=10, res=432, units='in')
        par(mfrow=c(2,2),mar=c(3,3,2,2),oma=c(3,3,2,2),pch=pch,cex=cex)
        plot(years,sapply(1:length(used_data), function(i) jack.res[[i]]$steepness$h),type="b",
             xlab="Year removed",ylab="",main="h in jackknife",ylim=resSR$steepness$h*ylim.range)
        abline(resSR$steepness$h,0,lwd=2,col=2,lty=2)

        plot(years,sapply(1:length(used_data), function(i) jack.res[[i]]$steepness$B0),type="b",
           xlab="Year removed",ylab="",main="B0 in jackknife",ylim=resSR$steepness$B0*ylim.range)
        abline(resSR$steepness$B0,0,lwd=2,col=2,lty=2)

        plot(years,sapply(1:length(used_data), function(i) jack.res[[i]]$steepness$R0),type="b",
         xlab="Year removed",ylab="",main="R0 in jackknife",ylim=resSR$steepness$R0*ylim.range)
        abline(resSR$steepness$R0,0,lwd=2,col=2,lty=2)

        plot(years,sapply(1:length(used_data), function(i) jack.res[[i]]$steepness$SB0),type="b",
           xlab="Year removed",ylab="",main="SB0 in jackknife",ylim=resSR$steepness$SB0*ylim.range)
        abline(resSR$steepness$SB0,0,lwd=2,col=2,lty=2)
        if (output) dev.off()
      }

    # plot SR curve ----
      if (output) png(file = paste0(filename,"_SRcurve.png"), width=8, height=5, res=432, units='in')
      par(mar=c(3,3,2,2),oma=c(3,3,2,2),mfrow=c(1,1))
      plot(data_SR$R ~ data_SR$SSB, cex=2, type = "p",xlab="SSB",ylab="R",pch=1,
           main="jackknife SR functions",ylim=c(0,max(data_SR$R)*1.3),xlim=c(0,max(data_SR$SSB)*1.3))
      points(rev(data_SR$SSB)[1],rev(data_SR$R)[1],col=1,type="p",lwd=3,pch=16,cex=2)
      for (i in 1:length(years)) {
        points(jack.res[[i]]$pred$SSB,jack.res[[i]]$pred$R,type="l",lwd=2,col=rgb(0,0,1,alpha=0.1))
      }
      points(resSR$pred$SSB,resSR$pred$R,col=2,type="l",lwd=3,lty=2)
      legend("topright", legend=c("Estimated SR function","A jackknife SR function"),col=c("red",rgb(0,0,1,alpha=0.1)),lty=c(2,1), cex=0.8)
      if (output) dev.off()
    }
    #fit.SRregime ----
    else {
      regime_unique = resSR$regime_resid$regime %>% unique()
      # plot parameters ----
      if (output) png(file = paste0(filename,"_pars.png"), width=15, height=7.5, res=432, units='in')
      par(mar=c(3,3,2,2),oma=c(3,3,2,2),mfrow=c(length(regime_unique),3),pch=pch,cex=cex)
      for(j in 1:length(regime_unique)) {
        plot(years,sapply(1:length(years), function(i) jack.res[[i]]$regime_pars$a[j]),type="b",
             xlab="Year removed",ylab="",
             main=paste0("a in jackknife (regime ",regime_unique[j],")"),ylim=resSR$regime_pars$a[j]*ylim.range)
        abline(resSR$regime_pars$a[j],0,lwd=2,col=2,lty=2)
        abline(v=resSR$input$regime.year-0.5,lwd=1,lty=3,col="blue")
        plot(years,sapply(1:length(years), function(i) jack.res[[i]]$regime_pars$b[j]),type="b",
             xlab="Year removed",ylab="",
             main=paste0("b in jackknife (regime ",regime_unique[j],")"),ylim=resSR$regime_pars$b[j]*ylim.range)
        abline(resSR$regime_pars$b[j],0,lwd=2,col=2,lty=2)
        abline(v=resSR$input$regime.year-0.5,lwd=1,lty=3,col="blue")
        plot(years,sapply(1:length(years), function(i) jack.res[[i]]$regime_pars$sd[j]),type="b",
             xlab="Year removed",ylab="",
             main=paste0("sd in jackknife (regime ",regime_unique[j],")"),ylim=resSR$regime_pars$sd[j]*ylim.range)
        abline(resSR$regime_pars$sd[j],0,lwd=2,col=2,lty=2)
        abline(v=resSR$input$regime.year-0.5,lwd=1,lty=3,col="blue")
      }
      if (output) dev.off()
      if(!is.null(resSR$steepness)){
        if (output) png(file = paste0(filename,"_steepness.png"), width=15, height=7.5, res=432, units='in')
        par(mar=c(3,3,2,2),oma=c(3,3,2,2),mfrow=c(length(regime_unique),4),pch=pch,cex=cex)
        for(j in 1:length(regime_unique)) {
        plot(years,sapply(1:length(years), function(i) jack.res[[i]]$steepness$h[j]),type="b",
             xlab="Year removed",ylab="",
             main=paste0("h in jackknife (regime ",regime_unique[j],")"),ylim=resSR$steepness$h[j]*ylim.range)
        abline(resSR$steepness$h[j],0,lwd=2,col=2,lty=2)
        abline(v=resSR$input$regime.year-0.5,lwd=1,lty=3,col="blue")

        plot(years,sapply(1:length(years), function(i) jack.res[[i]]$steepness$B0[j]),type="b",
             xlab="Year removed",ylab="",
             main=paste0("B0 in jackknife (regime ",regime_unique[j],")"),ylim=resSR$steepness$B0[j]*ylim.range)
        abline(resSR$steepness$B0[j],0,lwd=2,col=2,lty=2)
        abline(v=resSR$input$regime.year-0.5,lwd=1,lty=3,col="blue")
        plot(years,sapply(1:length(years), function(i) jack.res[[i]]$steepness$R0[j]),type="b",
             xlab="Year removed",ylab="",
             main=paste0("R0 in jackknife (regime ",regime_unique[j],")"),ylim=resSR$steepness$R0[j]*ylim.range)
        abline(resSR$steepness$R0[j],0,lwd=2,col=2,lty=2)
        abline(v=resSR$input$regime.year-0.5,lwd=1,lty=3,col="blue")
        plot(years,sapply(1:length(years), function(i) jack.res[[i]]$steepness$SB0[j]),type="b",
             xlab="Year removed",ylab="",
             main=paste0("SB0 in jackknife (regime ",regime_unique[j],")"),ylim=resSR$steepness$SB0[j]*ylim.range)
        abline(resSR$steepness$SB0[j],0,lwd=2,col=2,lty=2)
        abline(v=resSR$input$regime.year-0.5,lwd=1,lty=3,col="blue")
        }
        if (output) dev.off()
      }

      # plot SRcurve ----
      obs_data = resSR$pred_to_obs
      if (output) png(file = paste0(filename,"_SRcurve.png"), width=12, height=6, res=432, units='in')
      par(mar=c(3,3,2,2),oma=c(3,3,2,2),mfrow=c(1,length(regime_unique)))
      for(i in 1:length(regime_unique)) {
        use_data = dplyr::filter(obs_data,Regime == regime_unique[i])
        plot(use_data$R ~ use_data$SSB, cex=2, type = "p",pch=1,xlab="SSB",ylab="R",
             main=paste0("Jackknife SR for Regime ",regime_unique[i]),ylim=c(0,max(use_data$R)*1.3),xlim=c(0,max(obs_data$SSB)*1.3))
        for (j in 1:length(used_data)) {
          pred_data = jack.res[[j]]$pred %>% filter(Regime == regime_unique[i])
          points(pred_data$SSB,pred_data$R,type="l",lwd=2,col=rgb(0,0,1,alpha=0.1))
        }
        pred_data = resSR$pred %>% filter(Regime == regime_unique[i])
        points(pred_data$SSB,pred_data$R,col=2,type="l",lwd=3,lty=2)
       if(output) legend("topright", legend=c("Estimated SR function","A jackknife SR function"),col=c("red",rgb(0,0,1,alpha=0.1)),lty=c(2,1), cex=0.8)
       else if(i==1) legend("topleft", legend=c("Estimated SR function","A jackknife SR function"),col=c("red",rgb(0,0,1,alpha=0.1)),lty=c(2,1), cex=0.8)
      }
      if (output) dev.off()
    }
  }

  return(invisible(RES))
}

#' プロファイル尤度を計算する関数
#'
#' 結果のプロットもこの関数で行う
#' @param resSR \code{fit.SR}か\code{fit.SRregime}のオブジェクト
#' @param output pngファイルに出力するか否か
#' @param filename ファイル名
#' @param a_range 推定値に\code{a_range}を乗じた範囲で尤度計算を行う (\code{NULL}のとき\code{c(0.5,2)})
#' @param b_range 推定値に\code{a_range}を乗じた範囲で尤度計算を行う (\code{NULL}のとき\code{c(0.5,2)})
#' @param HS_b_restrict Hockey-StickのときにbをSSBの観測値の範囲にするか否か (\code{b_range}より優先される)
#' @param length 範囲を区切る数
#' @encoding UTF-8
#' @export
prof.likSR = function(resSR,output=FALSE,filename="Profile_Likelihood",a_range = NULL,b_range = NULL,HS_b_restrict = TRUE,length=50) {
  RES = list()
  if (is.null(a_range)) a_range = c(0.5,2)
  if (is.null(b_range)) b_range = c(0.5,2)
  if (class(resSR) == "fit.SR") {
    validate_sr(res_sr = resSR)
    a.grid <- seq(resSR$pars$a*a_range[1],resSR$pars$a*a_range[2],length=length)
    if (resSR$input$SR!="HS" || !isTRUE(HS_b_restrict)) {
      b.grid <- seq(resSR$pars$b*b_range[1],resSR$pars$b*b_range[2],length=length)
    } else {
      b.grid <- seq(min(resSR$input$SRdata$SSB),max(resSR$input$SRdata$SSB),length=length)
    }
    ba.grid = expand.grid(b=b.grid,a=a.grid)

    if (resSR$pars$rho==0 || resSR$input$out.AR==TRUE) {
      obj.f = function (a,b) resSR$obj.f(a=a,b=b,rho=0)
      prof.lik.res <- exp(-sapply(1:nrow(ba.grid), function(i) obj.f(a=ba.grid[i,2],b=ba.grid[i,1])))
    } else {
      obj.f = function (a,b,x) resSR$obj.f(a=a,b=b,rho=x)
      prof.lik.res <- exp(-sapply(1:nrow(ba.grid), function(i) {
        optimize(function(x) obj.f(a=ba.grid[i,2],b=ba.grid[i,1],x),interval=c(-0.999,0.999))$objective
      }))
    }
    if (output) png(file = paste0(filename,".png"), width=7.5, height=5, res=432, units='in')
    image(b.grid,a.grid,matrix(prof.lik.res,nrow=length),ann=F,col=cm.colors(12),
          ylim=range(a.grid),xlim=range(b.grid))
    par(new=T, xaxs="i",yaxs="i")
    contour(b.grid,a.grid,matrix(prof.lik.res,nrow=length),
            ylim=range(a.grid),xlim=range(b.grid),
            xlab="b",ylab="a",main="Profile Likelihood")
    points(resSR$pars$b,resSR$pars$a,cex=2,pch=4,col=2,lwd=3)
    if (output) dev.off()
    ba.grid.res = ba.grid
  } else { #fit.SRregime
    if (output) png(file = paste0(filename,".png"), width=15, height=5, res=432, units='in')
    par(mfrow=c(1,nrow(resSR$regime_pars)),mar=c(4,4,2,2))
    # par(mfrow=c(1,nrow(resSR$regime_pars)))
    prof.lik.res = NULL
    ba.grid.res = list()
    for (j in 1:nrow(resSR$regime_pars)) {
      a.grid <- seq(resSR$regime_pars$a[j]*a_range[1],resSR$regime_pars$a[j]*a_range[2],length=length)
      if (resSR$input$SR!="HS"|| !isTRUE(HS_b_restrict)) {
        b.grid <- seq(resSR$regime_pars$b[j]*b_range[1],resSR$regime_pars$b[j]*b_range[2],length=length)
      } else {
        if ("b" %in% resSR$input$regime.par) {
          ssb_range=range(resSR$pred_to_obs %>% dplyr::filter(Regime==resSR$regime_pars$regime[j]) %>%
                            dplyr::select(SSB))
          b.grid <- seq(ssb_range[1],ssb_range[2],length=length)
        } else {
          b.grid <- seq(min(resSR$input$SRdata$SSB),max(resSR$input$SRdata$SSB),length=length)
        }
      }
      ba.grid = expand.grid(b=b.grid,a=a.grid)
      ab_order = c("a"=which(resSR$regime_pars$a[j]==resSR$regime_pars$a),"b"=which(resSR$regime_pars$b[j]==resSR$regime_pars$b))
      a_fix = resSR$regime_pars$a[j]
      b_fix = resSR$regime_pars$b[j]
      ab = c(resSR$regime_pars$a,resSR$regime_pars$b)
      x = ab[!(ab %in% c(a_fix,b_fix))]
      obj.f = function(par_a,par_b,x) {
        a = resSR$regime_pars$a
        xa_length=length(a[-ab_order[1]])
        if (xa_length>0) a[-ab_order[1]] <- x[1:xa_length]
        b = resSR$regime_pars$b
        xb_length=length(b[-ab_order[2]])
        if (xb_length>0) b[-ab_order[2]] <- x[xa_length+(1:xb_length)]
        a[ab_order[1]] <- par_a
        b[ab_order[2]] <- par_b
        resSR$obj.f(a=a,b=b)
      }
      if (length(x)==1) {
        prof.lik.res <- cbind(prof.lik.res,exp(-sapply(1:nrow(ba.grid), function(i) {
          # opt = optim(x,obj.f,par_a=ba.grid[i,2],par_b=ba.grid[i,1],lower=x*1.0e-3,upper=x*1.0e+3,method="Brent")
          opt = optim(x,obj.f,par_a=ba.grid[i,2],par_b=ba.grid[i,1],method="BFGS")
          opt$value
        })))
      } else {
        prof.lik.res <- cbind(prof.lik.res,exp(-sapply(1:nrow(ba.grid), function(i) {
          # opt = optim(x,obj.f,par_a=ba.grid[i,2],par_b=ba.grid[i,1],lower=x*0.001,
          #             upper=x*1000,method="L-BFGS-B")
          opt = optim(x,obj.f,par_a=ba.grid[i,2],par_b=ba.grid[i,1],method="BFGS")
          opt$value
        })))
      }
      ba.grid.res[[j]] <- ba.grid

      image(b.grid,a.grid,matrix(prof.lik.res[,j],nrow=length),ann=F,col=cm.colors(12),
            ylim=range(a.grid),xlim=range(b.grid))
      par(new=T, xaxs="i",yaxs="i")
      contour(b.grid,a.grid,matrix(prof.lik.res[,j],nrow=length),
              ylim=range(a.grid),xlim=range(b.grid),
              xlab="b",ylab="a",main=paste0("Profile Likelihood for Regime ",resSR$regime_pars$regime[j]))
      points(resSR$regime_pars$b[j],resSR$regime_pars$a[j],cex=2,pch=4,col=2,lwd=3)
    }
    if (output) dev.off()
  }
  RES$prof.lik <- prof.lik.res
  RES$ba.grid <- ba.grid.res
  return(invisible(RES))
}

#' 再生産関係の推定結果をtxtファイルに出力する関数
#'
#' @param resSR \code{fit.SR}か\code{fit.SRregime}のオブジェクト
#' @param filename ファイル名('.txt')がつく
#' @encoding UTF-8
#' @export
out.SR = function(resSR,filename = "resSR") {
  RES = list()
  RES$SR = resSR$input$SR
  RES$method = resSR$input$method
  if (class(resSR) == "fit.SR") {
    validate_sr(res_sr = resSR)
    RES$AR = resSR$input$AR
    RES$out.AR = resSR$input$out.AR
    RES$pars = resSR$pars
  } else {
    RES$regime.year = resSR$input$regime.year
    RES$regime.key = resSR$input$regime.key
    RES$regime.par = resSR$input$regime.par
    RES$pars = resSR$regime_pars
  }
  RES$n = sum(resSR$input$w)
  RES$k = resSR$k
  RES$loglik = resSR$loglik
  RES$AIC = resSR$AIC
  if (!is.null(resSR$AIC.ar)) RES$AIC.ar = resSR$AIC.ar
  RES$AICc = resSR$AICc
  RES$BIC = resSR$BIC
  if(!is.null(resSR$input$bio_par)) RES$steepness=resSR$steepness
  RES$opt = resSR$opt
  if (class(resSR) == "fit.SR") {
    RES$pred_to_obs = as_tibble(resSR$input$SRdata) %>%
      dplyr::rename(Year = year) %>%
      dplyr::mutate(resid = resSR$resid,resid2 = resSR$resid2) %>%
      dplyr::mutate(Pred_from_SR = R/exp(resid),Pred_from_AR=R/exp(resid2)) %>%
      dplyr::select(Year,SSB,R,Pred_from_SR,resid,Pred_from_AR,resid2)
  } else {
    RES$pred_to_obs = resSR$pred_to_obs
  }
  RES$pred_to_obs =  as.data.frame(RES$pred_to_obs)

  capture.output(RES, file = paste0(filename,".txt"))
}

#' 再生産関係推定が収束しているかや最適解を得られているかを診断する関数
#'
#' @param resSR \code{fit.SR}か\code{fit.SRregime}のオブジェクト
#' @param n 初期値を変えてパラメータ推定する回数
#' @param sigma 初期値を変えるときの生起乱数の標準偏差
#' @param seed \code{set.seed}で使用するseed
#' @param output テキストファイルに結果を出力するか
#' @param filename ファイル名('.txt')がつく
#' @param fun_when_check5_replace check5にひっかかったとき、候補となるパラメータの中からどのパラメータを代表としてとってくるか。median, max, minなどの関数を入れる。デフォルト(今までの実装)はmedian。保守的な結果を出すならmaxのほうがよいかも(エクセルではmaxに近い値が出力されるみたい)
#' @encoding UTF-8
#' @export

check.SRfit = function(resSR,n=100,sigma=5,seed = 1,output=FALSE,filename="checkSRfit",
                       fun_when_check5_replace=median) {
  opt = resSR$opt
  SRdata = resSR$input$SRdata
  flag = rep(0,5)

  RES = list()
  # check convergence
  if (opt$convergence==0) {
    cat(RES$convergence <- "1. 収束しています (OK)","\n")
  } else {
    flag[1] <- 1
    stop(RES$convergence <- "1. 収束していないので初期値(引数p0)を変えて計算しなおしてください")
  }
  # hessian
  resSR2 = resSR
  if (is.null(resSR$opt$hessian)) {
    input = resSR$input
    input$p0 = resSR$opt$par
    input$hessian = TRUE
    if (class(resSR) == "fit.SR") {
      validate_sr(res_sr = resSR)
      #input$rep.opt = TRUE
      resSR2 = do.call(fit.SR, input)
    } else {
      if (class(resSR) == "fit.SRregime") {
        resSR2 = do.call(fit.SRregime, input)
      } else {
        input$rep.opt = TRUE
        resSR2 = do.call(fit.SRalpha, input)
      }
    }
  }
  if (all(diag(resSR2$opt$hessian) > 0)) {
    cat(RES$hessian <- "2. Hessian行列の対角成分が正定値になっています (OK)","\n")
  } else {
    flag[2] <- 1
    message(RES$hessian <- "2. Hessian行列の対角成分が正定値になっていません (HSかつ3以降のチェックがOKであれば問題ありません)")
  }

  # check boundary
  RES$boundary <- NULL
  if (class(resSR) == "fit.SR" || class(resSR) == "fit.SRalpha") {
    if (resSR$input$SR == "HS") {
      if (resSR$pars$b > max(SRdata$SSB)*0.99) message(RES$boundary <- "Parameter b reaching the maximum SSB")
      if (resSR$pars$b < min(SRdata$SSB)*1.01) message(RES$boundary <- "Parameter b reaching the minimum SSB")
    } else {
      if (1/resSR$pars$b > 10*max(SRdata$SSB)) message(RES$boundary <- "Proportional recruitment to SSB (no density-dependence)")
      if (1/resSR$pars$b < 0.1*min(SRdata$SSB)) message(RES$boundary <- "Extremely strong density-dependence of recruitment against SSB")
    }
  } else {
    for (i in 1:nrow(resSR$regime_pars)) {
      if ("b" %in% resSR$input$regime.par) {
        SRdata_r = dplyr::filter(resSR$pred_to_obs,Regime==resSR$regime_pars$regime[i])
      } else {
        SRdata_r = SRdata
      }
      if (resSR$input$SR == "HS") {
        if (resSR$regime_pars$b[i] > max(SRdata_r$SSB)*0.99) RES$boundary <- c(RES$boundary,paste0("Parameter b of regime ",resSR$regime_pars$regime[i], " reaching the maximum SSB"))
        if (resSR$regime_pars$b[i] < min(SRdata_r$SSB)*1.01) RES$boundary <- c(RES$boundary,paste0("Parameter b of regime ",resSR$regime_pars$regime[i], " reaching the minimum SSB"))
      } else {
        if (1/resSR$regime_pars$b[i] > 10*max(SRdata_r$SSB)) RES$boundary <- c(RES$boundary,paste0("Proportional recruitment of regime ",resSR$regime_pars$regime[i], " to SSB (no density-dependence)"))
        if (1/resSR$regime_pars$b[i] < 0.1*min(SRdata_r$SSB)) RES$boundary <- c(RES$boundary,paste0("Extremely strong density-dependence of recruitment against SSB in ",resSR$regime_pars$regime[i]))
      }
    }
  }
  if (is.null(RES$boundary)) {
    cat(RES$boundary <- "3. どの推定パラメータも壁(boundaries)にあたっていないのでOKです (OK)","\n")
  } else {
    flag[3] <- 1
    for (i in 1:length(RES$boundary)) {
      message("3. パラメータは壁にあたっています(HSで折れ点が過去最小・最大親魚量になっているときにそうなります。HSでない場合は推定の不安定性を示唆します。)", RES$boundary[i])
    }
  }

  set.seed(seed)
  loglik = NULL
  pars = NULL
  resSR_list = list()
  for (i in 1:n) {
    input = resSR$input
    for (j in 1:100) {
      input$p0 <- opt$par + rnorm(length(opt$par),0,sigma)
      if (class(resSR) == "fit.SR") {
        #input$rep.opt = TRUE
        resSR2 = try(do.call(fit.SR, input),silent=TRUE)
      } else {
        if (class(resSR) == "fit.SRregime") {
          resSR2 = try(do.call(fit.SRregime, input),silent=TRUE)
        } else {
          input$rep.opt = TRUE
          resSR2 = try(do.call(fit.SRalpha, input),silent=TRUE)
        }
      }
      if (class(resSR2) != "try-error") break
    }
    resSR_list[[i]] = resSR2
    loglik = c(loglik, resSR2$loglik)
    pars = rbind(pars,resSR2$opt$par)
  }
  max_loglik = max(loglik)
  optimal = NULL
  if (resSR$loglik-max_loglik < -0.001) {
    flag[4] <- 1
    diff_loglik = abs(resSR$loglik-max_loglik)
    message(paste0("Maximum difference of log-likelihood is ",round(diff_loglik,6)))
    optimal = resSR_list[[which.max(loglik)]]
    message(RES$optim <- str_c("4. パラメータが大域解に達していません (fit.SRtolで再計算をおこなうか、本関数の返り値のoptimumに結果を置き換えてください。大域解を得るための初期値は,",str_c(optimal$input$p0,collapse=", "),"です)"))
    RES$loglik_diff = diff_loglik
  } else {
    cat(RES$optim <- "4. パラメータが大域解に達しているのでOKです (OK)","\n")
    # global optimumに達している場合のみ
    loglik_diff = purrr::map_dbl(loglik, function(x) abs(diff(c(x,max_loglik))))
    problem = NULL
    diff_threshold = 1.0e-6
    # a_diff = NULL; b_diff = NULL; sd_diff = NULL; rho_diff = NULL
    for (i in 1:n) {
      if (loglik_diff[i] < diff_threshold) {
        if (all(abs(pars[i,] - resSR$opt$par) < 0.001)) {
          problem = c(problem,FALSE)
          # a_diff = c(a_diff,0); b_diff = c(b_diff,0); sd_diff = c(sd_diff,0)
          # if (class(resSR)=="fit.SR" && resSR$pars$rho != 0) rho_diff = c(rho_diff,NULL)
        } else {
          problem = c(problem,TRUE)
          # a_diff = c(a_diff,max(abs(resSR_list[[i]]$pars$a/resSR$pars$a-1))*100)
          # b_diff = c(b_diff,max(abs(resSR_list[[i]]$pars$b/resSR$pars$b-1))*100)
          # sd_diff = c(sd_diff,max(abs(resSR_list[[i]]$pars$sd/resSR$pars$sd-1))*100)
          # if (class(resSR)=="fit.SR" && resSR$pars$rho != 0) {
          #   rho_diff = c(rho_diff,max(abs(resSR_list[[i]]$pars$rho/resSR$pars$rho-1))*100)
          # }
        }
      } else {
        problem = c(problem,FALSE)
      }
    }
    if (sum(problem)>0) {
      flag[5] <- 1
      RES$loglik_diff <- loglik_diff
      message(RES$pars <- str_c("5. 同じ最大尤度(",diff_threshold,"よりも小さい違い)を持つ複数のパラメータが見つかりました(L1かつHSでよく見られます)。"))
      # RES$percent_bias = c("a"=max(a_diff),"b"=max(b_diff),"sd" = max(sd_diff))
      # message("Maximum percent bias of 'a' is ", round(as.numeric(RES$percent_bias["a"]),6),"%")
      # message("Maximum percent bias of 'b' is ", round(as.numeric(RES$percent_bias["b"]),6),"%")
      # message("Maximum percent bias of 'sd' is ", round(as.numeric(RES$percent_bias["sd"]),6),"%")
      # if (class(resSR)=="fit.SR" && resSR$pars$rho != 0) {
      #   RES$percent_bias = c(RES$percent_bias,"rho" = max(rho_diff))
      #   message("Maximum percent bias of 'rho' is ", round(as.numeric(RES$percent_bias["rho"]),6),"%")
      # }
      selected_set <- which(loglik_diff<diff_threshold)
      if(class(resSR)!="fit.SRregime"){
        par_list = t(sapply(1:n, function(i) unlist(resSR_list[[i]]$pars)[unlist(resSR$pars) != 0]))
        par_list0 <- par_list
        par_list = par_list[selected_set,]
        bias_list = t(sapply(1:n, function(i) 100*(unlist(resSR_list[[i]]$pars)[unlist(resSR$pars) != 0]/unlist(resSR$pars)[unlist(resSR$pars)!=0]-1)))
        bias_list = bias_list[selected_set,]
      }else{
        par_list = t(sapply(1:n, function(i) unlist(resSR_list[[i]]$regime_pars)[unlist(resSR$regime_pars) != 0]))
        par_list0 <- par_list
        par_list = par_list[selected_set,]
        bias_list = t(sapply(1:n, function(i) 100*(unlist(resSR_list[[i]]$regime_pars)[unlist(resSR$regime_pars) != 0]/unlist(resSR$regime_pars)[unlist(resSR$regime_pars)!=0]-1)))
        bias_list = bias_list[selected_set,]
      }
      par_summary = apply(par_list,2,summary)
      percent_bias_summary = apply(bias_list,2,summary)
      RES$par_summary <- par_summary
      RES$percent_bias_summary <- percent_bias_summary
      RES$par_list <- par_list
      RES$par_list0 <- par_list0
      # すべてのパラメータのfun_when_check5_replaceに最も近いパラメータセットを持つものを選んでoptimalに入れちゃう
      #      x <- sweep(par_list,2,apply(par_list,2,median),FUN="/") %>% apply(1,mean)
     if(class(resSR)!="fit.SRregime"){
        x <- sweep(par_list[,c("a","b")],2,apply(par_list[,c("a","b")],2,fun_when_check5_replace),FUN="/") %>% apply(1,mean)
     }else{
        tmp <- 2:(1+2*length(unique(resSR$input$regime.key)))
        x <- sweep(par_list[,tmp],2,apply(par_list[,tmp],2,fun_when_check5_replace),FUN="/") %>% apply(1,mean)
      }
      selected <- which.min(abs(x-1))
      cat("ほとんど同じ尤度を持つパラメータの範囲 (",n,"回試行のうち",nrow(par_list),"回分),\n")
      print(apply(par_list,2,summary))
      optimal <- resSR_list[selected_set][[selected]]
      cat(round(as.numeric(unlist(optimal$pars)),4), "をoptimumに出力します(そのときの初期値は",str_c(optimal$input$p0,collapse="-"),"です)\n")
    } else {
      cat(RES$pars <- "5. パラメータが唯一の解として推定されているのでOKです (OK)","\n")
    }
  }
  if (output) {
    capture.output(RES,file=paste0(filename,".txt"))
  }
  if (!is.null(optimal)) RES$optimum = optimal
  # RES$loglik = loglik
  # RES$par_list = par_list
  # RES$percent_bias_list = bias_list
  RES$flag <- flag
  RES$loglik <- loglik
  return(RES)
}

#' 再生産関係推定でSR=HSにおいて同一尤度を持つ複数のパラメータセットがえられたとき、SR=Mesnilとしてcheck.SRfitを呼び出し、gammaを変えながら尤度が一意に定まる値を探す関数
#'
#' @param resSR \code{fit.SR}か\code{fit.SRregime}のオブジェクトで、かつSR=Mesnilのもの
#' @param n 初期値を変えてパラメータ推定する回数 (check.SRfitの引数)
#' @param sigma 初期値を変えるときの生起乱数の標準偏差(check.SRfitの引数)
#' @param seed \code{set.seed}で使用するseed(check.SRfitの引数)
#' @param gamma_ini SR=Mesnilの曲がり具合を決めるパラメータ デフォルトは10で、9,8,…,2,1,0.9,0.8,…と複数のパラメータが得られるまで探していく。
#' @encoding UTF-8
#' @export
specify.Mesnil.gamma <- function(resSR,n=100,seed = 1,sigma=5,gamma_ini=10){
  opt = resSR$opt
  SRdata = resSR$input$SRdata
  resSR2 = resSR
  RES = list()

  if (resSR$input$SR!="Mesnil") stop("This function checks convergence of the fitting SR function if SR=Mesnil.")

  identical.likely<-TRUE
  gamma_post <- gamma_ini
  while(identical.likely){
    set.seed(seed)
    loglik = NULL
    pars = NULL
    resSR_list = list()
    for (i in 1:n) {
      input = resSR$input
      input$gamma = gamma_post

      for (j in 1:100) {
        input$p0 <- opt$par + rnorm(length(opt$par),0,sigma)
        if (class(resSR) == "fit.SR") {
          #input$rep.opt = TRUE
          resSR2 = try(do.call(fit.SR, input),silent=TRUE)
        } else {
          if (class(resSR) == "fit.SRregime") {
            resSR2 = try(do.call(fit.SRregime, input),silent=TRUE)
          } else {
            input$rep.opt = TRUE
            resSR2 = try(do.call(fit.SRalpha, input),silent=TRUE)
          }
        }
        if (class(resSR2) != "try-error") break
      }
      resSR_list[[i]] = resSR2
      loglik = c(loglik, resSR2$loglik)
      pars = rbind(pars,resSR2$opt$par)
    }
    max_loglik = max(loglik)
    optimal = NULL
    if (resSR$loglik-max_loglik < -0.001) {
      diff_loglik = abs(resSR$loglik-max_loglik)
      message(paste0("Maximum difference of log-likelihood is ",round(diff_loglik,6)))
      optimal = resSR_list[[which.max(loglik)]]
      RES$loglik_diff = diff_loglik
    } else {
      # global optimumに達している場合のみ
      loglik_diff = purrr::map_dbl(loglik, function(x) abs(diff(c(x,max_loglik))))
      problem = NULL
      diff_threshold = 1.0e-6
      # a_diff = NULL; b_diff = NULL; sd_diff = NULL; rho_diff = NULL
      for (i in 1:n) {
        if (loglik_diff[i] < diff_threshold) {
          if (all(abs(pars[i,] - resSR$opt$par) < 0.001)) {
            problem = c(problem,FALSE)
          } else {
            problem = c(problem,TRUE)
          }
        } else {
          problem = c(problem,FALSE)
        }
      }
      if (sum(problem)>0) {
        RES$loglik_diff <- loglik_diff
        message(RES$pars <- str_c("同じ最大尤度(",diff_threshold,"よりも小さい違い)を持つ複数のパラメータが見つかりました(L1かつHSでよく見られます)。"))
        selected_set <- which(loglik_diff<diff_threshold)
        if(class(resSR)!="fit.SRregime"){
          par_list = t(sapply(1:n, function(i) unlist(resSR_list[[i]]$pars)[unlist(resSR$pars) != 0]))
          par_list0 <- par_list
          par_list = par_list[selected_set,]
          bias_list = t(sapply(1:n, function(i) 100*(unlist(resSR_list[[i]]$pars)[unlist(resSR$pars) != 0]/unlist(resSR$pars)[unlist(resSR$pars)!=0]-1)))
          bias_list = bias_list[selected_set,]
        }else{
          par_list = t(sapply(1:n, function(i) unlist(resSR_list[[i]]$regime_pars)[unlist(resSR$regime_pars) != 0]))
          par_list0 <- par_list
          par_list = par_list[selected_set,]
          bias_list = t(sapply(1:n, function(i) 100*(unlist(resSR_list[[i]]$regime_pars)[unlist(resSR$regime_pars) != 0]/unlist(resSR$regime_pars)[unlist(resSR$regime_pars)!=0]-1)))
          bias_list = bias_list[selected_set,]
        }
        par_summary = apply(par_list,2,summary)
        percent_bias_summary = apply(bias_list,2,summary)
        RES$par_summary <- par_summary
        RES$percent_bias_summary <- percent_bias_summary
        RES$par_list <- par_list
        RES$par_list0 <- par_list0
        # すべてのパラメータのmedianに最も近いパラメータセットを持つものを選んでoptimalに入れちゃう
        #      x <- sweep(par_list,2,apply(par_list,2,median),FUN="/") %>% apply(1,mean)
        if(class(resSR)!="fit.SRregime"){
          x <- sweep(par_list[,c("a","b")],2,apply(par_list[,c("a","b")],2,median),FUN="/") %>% apply(1,mean)
        }else{
          tmp <- 2:(1+2*length(unique(resSR$input$regime.key)))
          x <- sweep(par_list[,tmp],2,apply(par_list[,tmp],2,median),FUN="/") %>% apply(1,mean)
        }
        selected <- which.min(abs(x-1))
        optimal <- resSR_list[selected_set][[selected]]

        identical.likely<-FALSE
      }else{
        gamma_pre <- input$gamma
        degit_n <- floor(log10(gamma_pre))
        ssd <- gamma_pre %/% 10^(degit_n)
        if(ssd!=1){
          gamma_post <- gamma_pre - 10^(degit_n)
        }else{
          gamma_post <- gamma_pre - 10^(degit_n-1)
        }
        input$gamma <- gamma_post
      }
    }

  }

  if(gamma_pre==gamma_ini)stop("Mesnil with gamma_ini has multi-identical likelyhood. Set gamma_ini larger than current value.")
  input$gamma = gamma_pre
  if (class(resSR) == "fit.SR") {
    #input$rep.opt = TRUE
    optimal = try(do.call(fit.SR, input),silent=TRUE)
  } else {
    if (class(resSR) == "fit.SRregime") {
      optimal = try(do.call(fit.SRregime, input),silent=TRUE)
    } else {
      input$rep.opt = TRUE
      optimal = try(do.call(fit.SRalpha, input),silent=TRUE)
    }
  }
  print(paste0("unique SR parameter was given where gamma = ",gamma_pre,"."))
  RES$optimum = optimal
  RES$loglik <- loglik
  RES$gamma <- gamma_pre
  return(RES)
}


#' 再生産関係のブートストラップの結果をggplotで生成する関数する関数
#'
#' @param boot.res \code{boot.SR}のオブジェクト
#' @param CI プロットする信頼区間
#' @encoding UTF-8
#' @export

bootSR.ggplot = function(boot.res, CI=0.80) {
  if (class(boot.res$input$Res) == "fit.SRregime") {
    stop("This function cannot handle 'fit.SRregime' at present")
  }
  estimate_data = boot.res$input$Res$pred
  obs_data = as_tibble(boot.res$input$Res$input$SRdata)

  for (i in 1:boot.res$input$n) {
    tmp_tbl = boot.res[[i]]$pred %>% mutate(simID = i)
    if (i == 1) {
      boot_pred = tmp_tbl
    } else {
      boot_pred = bind_rows(boot_pred,tmp_tbl)
    }
  }

  summary_boot = boot_pred %>%
    group_by(SSB) %>%
    summarise(median = median(R), lower=quantile(R, probs=0.5*(1-CI))[1], upper=quantile(R, probs=1-0.5*(1-CI))[1])

  g1 = ggplot(data = NULL)+
    geom_ribbon(data = summary_boot, aes(x=SSB,ymin=lower,ymax=upper),alpha=0.3,fill="red")+
    geom_path(data = estimate_data,aes(x=SSB,y=R),size=2) +
    geom_point(data = obs_data,aes(x=SSB,y=R),size=2)+
    geom_path(data = summary_boot, aes(x=SSB,y=median),colour="red",size=2,linetype="dashed")+
    theme_SH()
  g1
}

#' 再生産関係の推定パラメータの相関を出力する関数
#'
#' @inheritParams fit.SR
#' @inheritParams fit.SRregime
#' @param resSR \code{fit.SR}または\code{fit.SRregime}のオブジェクト
#' @return 以下の要素からなるリスト
#' \describe{
#' \item{\code{hessian}}{ヘッセ行列}
#' \item{\code{cov}}{推定されたパラメータの分散共分散行列}
#' \item{\code{cor}}{推定されたパラメータの相関行列}
#' }
#' @examples
#' \dontrun{
#' data(res_vpa)
#' SRdata <- get.SRdata(res_vpa)
#' resSR <- fit.SR(SRdata, SR = c("HS","BH","RI")[1],
#'                 method = c("L1","L2")[2], AR = 1,
#'                 out.AR = TRUE)
#' corRes = corSR(resSR)
#' corRes$cor
#' }
#' @encoding UTF-8
#' @export
corSR = function(resSR) {
  if (!resSR$input$hessian) {
    resSR$input$hessian <- TRUE
    resSR$input$p0 = resSR$opt$par
    if (class(resSR) == "fit.SR") resSR = do.call(fit.SR, resSR$input)
    if (class(resSR) == "fit.SRregime") resSR = do.call(fit.SRregime, resSR$input)
  }
  hessian = resSR$opt$hessian
  cov = solve(hessian)
  cor = stats::cov2cor(cov)
  return (list(hessian=hessian,cov=cov,cor=cor))
}

#' Steepness (h) と関連するパラメータ (SB0,R0,B0)を計算する関数
#'
#' @param SR "HS", "BH", "RI"のいずれか
#' @param rec_pars 再生産関係のパラメータで\code{rec_pars$a},\code{rec_pars$b}で使用する
#' @param M 年齢別自然死亡係数 (ベクトルで与えるか、年齢共通の場合\code{M=0.4}のようにしてもよい)
#' @param waa (親魚量の)年齢別体重
#' @param maa 年齢別親魚量
#' @param plus_group 最高齢がプラスグループかどうか
#' @return 以下の要素からなるデータフレーム
#' \describe{
#' \item{\code{SPR0}}{F=0のときのSPR(この逆数がreplacement lineの傾き)}
#' \item{\code{SB0}}{F=0のときの親魚量}
#' \item{\code{R0}}{F=0のときの加入量}
#' \item{\code{B0}}{F=0のときの資源量}
#' \item{\code{h}}{steepness: BHかRIのときは0.2×SB0のときの加入量がh×R0, HSのときはh=1-b/SB0}
#' }
#' @examples
#' \dontrun{
#' data(res_vpa)
#' SRdata <- get.SRdata(res_vpa)
#' resSR <- fit.SR(SRdata, SR = c("HS","BH","RI")[1],
#'                 method = c("L1","L2")[2], AR = 1,
#'                out.AR = TRUE)
#' rec_pars = resSR$pars
#' year <- "2017"
#' M = res_vpa$input$dat$M[,year]
#' waa = res_vpa$input$dat$waa[,year]
#' maa = res_vpa$input$dat$maa[,year]
#' Res_h = calc_steepness(SR="HS",rec_pars=rec_pars,M=M,waa=waa,maa=maa,plus_group=TRUE)
#' Res_h
#' }
#' @encoding UTF-8
#' @export
calc_steepness = function(SR="HS",rec_pars,M,waa,maa,plus_group=TRUE,faa = NULL, Pope=TRUE) {
  if (length(M)==1) {
    M = rep(M,length(waa))
  }
  if (length(waa) != length(maa) || length(M) != length(maa)) {
    stop("The lengths of 'waa' and 'maa' must be equal")
  }
  is_MSY <- ifelse(is.null(faa),0,1)
  if (is.null(faa)) {
    faa <- rep(0,length(waa))
  }
  est_detMSY = function(x) {
    NAA0 = 1
    for (i in 1:(length(waa)-1)) {
      NAA0 = c(NAA0,rev(NAA0)[1]*exp(-M[i]-x*faa[i]))
    }
    if (plus_group) {
      NAA0[length(NAA0)] = rev(NAA0)[1]/(1-exp(-1*rev(M)[1]-x*rev(faa)[1]))
    }
    BAA0 = NAA0*waa
    SSB0 = BAA0*maa
    SPR0 = sum(SSB0) #get.SRRと一致 (testに使える)

    # 再生産関係とy=(1/SPR0)*xの交点を求める
    rec_a = rec_pars$a
    rec_b = rec_pars$b
    validate_sr(SR = SR)
    if (SR == "HS") {
      R0 = rec_pars$a*rec_pars$b
      SB0 = R0*SPR0
      if (SB0<rec_b) {
        # warning("Virgin equilibrium does not exist!")
        R0 <- 0
        SB0 <- 0
      }
      h = (SB0-rec_b)/SB0
    }
    if (SR == "Mesnil") {
      gamma <- rec_pars$gamma
      K = sqrt(rec_b^2+gamma^2/4)
      SB0 = (2*K/(SPR0*rec_a/2)-2*rec_b-2*K)/(1/(SPR0*rec_a/2)^2-2/(SPR0*rec_a/2)) #Mesnil and Rochet 2010 ICES JMSより
      R0 = SB0/SPR0
      if (1/SPR0>rec_a) {
        # warning("Virgin equilibrium does not exist!")
        R0 <- 0
        SB0 <- 0
      }
      h = (SB0-rec_b)/SB0 # steepnessの定義はHSと同じでよい?
    }
    if (SR == "BH") {
      SB0 = (rec_a*SPR0-1)/rec_b
      R0 = SB0/SPR0
      h = (rec_a*0.2*SB0/(1+rec_b*0.2*SB0))/R0
    }
    if (SR == "RI") {
      SB0 = (1/rec_b)*log(rec_a*SPR0)
      R0 = SB0/SPR0
      h = (rec_a*0.2*SB0*exp(-rec_b*0.2*SB0))/R0
    }
    B0 = sum(R0*BAA0)
    Res = data.frame(SPR0 = SPR0, SB0 = SB0, R0 = R0, B0 = B0, h = h)

    if(is_MSY==1){
      ypr.spr = ref.F(Fcurrent=x*faa,M=M,waa=waa,waa.catch = waa,maa =maa,
                      Pope=Pope,pSPR=NULL,F.range=NULL,plot=FALSE)
      ypr.spr <- ypr.spr$ypr.spr[1,]
      Yield <- as.numeric(R0*ypr.spr["ypr"])
      Res = cbind(Res,data.frame(Yield=Yield))
    }
    return(Res)
  }

  RES = est_detMSY(0)
  if (is_MSY==1) RES = RES %>% dplyr::select(-Yield)
  if(is_MSY==1){
    obj_fun = function(x) {
      RES = est_detMSY(x)
      return(-RES$Yield)
    }
    x_grid = seq(0,10,length=101)
    tmp = x_grid %>% purrr::map_dbl(.,obj_fun)
    Opt = optimize(obj_fun,x_grid[c(max(1,which.min(tmp)-1),min(which.min(tmp)+1,length(x_grid)))])
    Fmsy2F=Opt$minimum
    RES2 = est_detMSY(Fmsy2F)
    RES2 = RES2 %>% dplyr::select(-h) %>% dplyr::mutate(Fmsy2F=Fmsy2F)
    colnames(RES2) <- c("SPRmsy","SBmsy","Rmsy","Bmsy","MSY","Fmsy2F")
    RES = cbind(RES,RES2)
  }
  return(RES)
}


#' 再生産関係の推定結果からbootstrapを実行し、hの分布を計算する関数
#'
#' @param res_SR fit.SRまたはfit.SRregime、モデル平均の結果
#' @param M 年齢別自然死亡係数 (ベクトルで与えるか、年齢共通の場合\code{M=0.4}のようにしてもよい)
#' @param waa (親魚量の)年齢別体重
#' @param maa 年齢別親魚量
#' @param plus_group 最高齢がプラスグループかどうか
#'
#' @export
#'


boot_steepness <- function(res_SR, M, waa, maa, n=100, plus_group=TRUE){

    is.model.average <- class(res_SR)!="fit.SRregime" && class(res_SR)!="fit.SR"

    if(is.model.average){ # model average case (calculate recursively)
      res_steepness <- purrr::map_dfr(res_SR, boot_steepness, n=n, M=M, waa=waa, maa=maa, plus_group=plus_group, .id="id")
    }else{
      res_boot <- boot.SR(res_SR, n=n)
      if(class(res_boot[[1]])=="fit.SRregime"){ # regime shift
          res_steepness <- purrr::map_dfr(res_boot[1:n], function(x){
              par.matrix <- x$regime_pars[c("a","b")]
              tmplist <- purrr::map_dfr(seq_len(nrow(par.matrix)),
                                    function(i){
                                        calc_steepness(SR=res_SR$input$SR,rec_pars=par.matrix[i,],M=M,waa=waa,maa=maa,
                                                       plus_group=plus_group)
                                    },.id="id")
              tmplist <- bind_cols(tmplist,par.matrix)
          })
      }
      if(class(res_boot[[1]])=="fit.SR"){ # normal
          res_steepness <- purrr::map_dfr(res_boot[1:n], function(x){
              calc_steepness(SR=res_SR$input$SR,rec_pars=x$pars,M=M,waa=waa,maa=maa,plus_group=plus_group)
          })
          res_steepness$id <- 1
          res_steepness <- res_steepness %>%
              bind_cols(purrr::map_dfr(res_boot[1:n], function(x) x$pars))
      }}

    res_steepness
}

#' fit.SRとfit.SRregimeの重み付け(どのデータを使うか)の設定方法は、引数wで与える場合とSRdata$weightで与える場合の2パターンある。両者が矛盾なく設定されているかを確認するための関数

check_consistent_w <- function(w, SRdata){
    if(is.null(SRdata$weight)  && is.null(w)) SRdata$weight <- w <- rep(1,length(SRdata$R))
    if(is.null(SRdata$weight)  && !is.null(w)) SRdata$weight <- w
    if(!is.null(SRdata$weight) &&  is.null(w)) w <- SRdata$weight
    if(!is.null(SRdata$weight) && !is.null(w) && sum(unlist(SRdata$weight)!=unlist(w))>0) stop("SRdata$weight と引数wの両方に重み付けの指定がなされていて、両者が違います")
    lst(w,SRdata)
}


#'
#' 再生産関係を網羅的にフィットする
#'
#' regimeなしのものだけに対応
#'
#' @export
#'

tryall_SR <- function(data_SR, plus_group=TRUE, bio_par=NULL, tol=FALSE){

  SRmodel.list <- expand.grid(SR.rel = c("HS","BH","RI"), L.type = c("L1", "L2")) %>%
    as_tibble()

  if(tol==TRUE) fit.SR <- fit.SR_tol
  SRmodel.list$pars <- purrr::map2(SRmodel.list$SR.rel, SRmodel.list$L.type,
                                   function(x,y){
                                     res1 <- unlist(fit.SR(data_SR, SR = x, method = y,
                                                           AR = 0, hessian = FALSE, out.AR=TRUE,
                                                           bio_par = bio_par, plus_group = plus_group)
                                                    [c("pars","AICc","steepness")])
                                     tmp <- fit.SR(data_SR, SR = x, method = y,
                                                   AR = 1, hessian = FALSE, out.AR=TRUE,
                                                   bio_par = bio_par, plus_group = plus_group)
                                     res2 <- unlist(tmp[c("pars","AICc","steepness")])
                                     res2 <- c(res2,"deltaAIC(AIC_AR-AIC_noAR)"=tmp$AIC.ar[2]-tmp$AIC.ar[1])
                                     res3 <- unlist(fit.SR(data_SR, SR = x, method = y,
                                                           AR = 1, hessian = FALSE, out.AR=FALSE,
                                                           bio_par = bio_par, plus_group = plus_group)[c("pars","AICc","steepness")])
                                     bind_rows(res1,res2,res3,.id="id")
                                   })

  SRmodel.list <- SRmodel.list %>%
    unnest(col=pars) %>%
    left_join(tibble(id=as.character(1:3),AR.type=c("non","outer","inner"))) %>%
    arrange(AICc,AR.type)

  return(SRmodel.list)

}

#'
#' 再生産関係のフィット(頑健版)
#'
#' fit.SRのwrapper。与えられたinputはそのままfit.SRに渡すが、計算したあとにcheck.SRfitを回し、大域解に達していない場合や同じ尤度を持つ複数の大域解がある場合、その中央値をとるres_SRを返す
#'
#' @param ... fit.SRへの引数
#' @param n_check check.SRfitに渡す引数n(初期値を変えて計算を繰り返す回数)
#' @param seed check.SRfitに渡す引数seed。初期値を変えるときに使う乱数。toolのスクリプトではずっと12345を利用していた関係で、デフォルトの引数は12345にしている。
#'
#' @details fit.SR_tolの動作は以下の通りです。1) fit.SRを実施 2) fit.SRの結果(便宜的にres0とする)にcheck.SRfitをあてはめる。check.SRfitでは以下のことが実施されます。2-1) 初期値を変えてn回パラメータ推定を繰り返す。2-2) res0から得られている対数尤度と、n回パラメータ推定を繰り返したときに得られた対数尤度を比較し、res0の対数尤度よりも大きいものが見つかったらres0は「4. 大域解ではない」という判定をする。また、帰り値$optimumの中に、最大の対数尤度を得たときの推定結果を返す。2-3) 4.がOKの場合、n回分の対数尤度の中で最大値との差が1-06以下だけど、パラメータの推定値が異なる(0.001以上)場合には、「5. 同じ対数尤度を持つ複数のパラメータが見つかりました」というメッセージが出る。この場合、対数尤度の差が1-06以下を示すパラメータの範囲(最小、最大、平均など)が出力される。また、中央値に最も近いパラメータのセットを$optimumに返す。3) check.SRfitの返り値を見て、4のフラグが立っている場合、res0をcheck.SRfitの返り値$optimumに置き換える。そうなった場合、もう一度check.SRfitを実行する(これによって4のフラグは立たなくなることが想定)。4) check.SRfitの返り値を見て、5のフラグが立っている場合、res0をcheck.SRfitの返り値$optimumに置き換える。5) res0を返す。
#'
#' @export
#'

fit.SR_tol <- function(...,n_check=100,is_regime=FALSE,seed=12345,fun_when_check5_replace=median){

  if(is_regime==FALSE)  res_SR <- fit.SR(...)
  if(is_regime==TRUE )  res_SR <- fit.SRregime(...)

  cat("...check.SRfit 1 回目....\n")
  check <- check.SRfit(res_SR, output=FALSE, n=n_check, fun_when_check5_replace=fun_when_check5_replace)
  # まず大域解に達しているかどうかのflagを確認し、flagが立っている場合には$optimと置き換える
  if(check$flag[4]==1) {
    res_SR <- check$optimum
    cat("大域解を得るための初期値に変えたres_SRの結果に置き換えます\n")

    # 再度check.SRfit
    cat("...check.SRfit 2回目....\n")
    check <- check.SRfit(res_SR,output=FALSE,n=n_check, fun_when_check5_replace=fun_when_check5_replace)
  }

  # 5番目のflag(尤度が同じパラメータの範囲が存在する)が立っているかどうかも確認し、flagが立っていればoptimumと置き換える
  if(check$flag[5]==1){
    res_SR <- check$optimum
    cat("尤度が同じパラメータの範囲で中央値をとる結果に置き換えます\n")
  }
  return(res_SR)
}


#'
#' fit.SRのwrapper関数でhに制約を置いたもの
#'
#' 現状は推定されたパラメータのみが出力される
#'
#' @export
#'
#'

fit.SR_pen <- function(bio_par, h_upper=Inf, h_lower=0.2, plus_group=TRUE, ...){
    res1 = fit.SR(...)
#    ri1$pars
#    ri1$steepness
#    h_upper = 1
    init = res1$opt$par

    obj_pen = function(x,out=FALSE) {
        a = exp(x[1]); b= exp(x[2]);
        h0 = calc_steepness(SR="RI",rec_pars=data.frame("a"=a,"b"=b),waa=bio_par$waa,maa=bio_par$maa,M=bio_par$M,plus_group=plus_group)
        h <- h0[1,"h"]
        if(out==FALSE) {
            res1$obj.f2(x)+1000*max(0,h-h_upper) + 1000*max(0,h_lower-h)
        } else{
            bind_cols(data.frame("a"=a,"b"=b),h0)
        }
    }

    opt = optim(init,obj_pen)
    obj_pen(opt$par,out=TRUE)

}


#'
#' 再生産関係のレジームシフトを隠れマルコフモデルで推定する
#'
#' Tang et al. (2021) ICESJMSの論文を改変
#'
#' @inheritParams fit.SRregime
#' @params k_regime 推定するレジームの数(2以上の整数)
#' @params b_range パラメータbの範囲
#' @params p0 パラメータの初期値(リスト)
#' @params overwrite cppファイルのコンパイルを上書きして行うか
#' @encoding UTF-8
#' @export
#'
#
hmm_SR = function(SRdata,SR="BH",k_regime=2,gamma=0.01,b_range=NULL,p0=NULL,overwrite=FALSE,max.ssb.pred = 1.3) {
  argname <- ls()
  arglist <- lapply(argname,function(xx) eval(parse(text=xx)))
  names(arglist) <- argname

  k_regime<-round(k_regime)

  if(k_regime<2) stop("Incorrect 'k_regime'")

  data = SRdata
  st = data$SSB
  # st = st/(10^4)
  yt = log(data$R/data$SSB)

  SRcode = case_when(SR=="RI" ~ 1,SR=="BH" ~ 2,SR=="HS"~3,SR=="Mesnil"~4,TRUE~5)
  if (SRcode==5) stop("SR not recognized")
  if (is.null(b_range)) {
    if (SRcode<3) { # Ricker or BH
      b_range = range(1/st)
    } else {
      b_range = range(st)
    }
  }

  tmb_data = list(st=st,yt=yt,
                  alpha_u=max(yt),alpha_l=min(yt),
                  beta_u=max(b_range),beta_l=min(b_range),
                  sigma_u=sd(yt),SRcode=SRcode,gamma=gamma)

  if (!is.null(p0)) {
    parameters = p0
  } else {
    parameters = list(
      lalpha = -log(k_regime+1-1:k_regime),
      lbeta = rep(0,k_regime),
      lsigma = rep(0,k_regime),
      pi1_tran = rep(0,k_regime-1),
      qij_tran = matrix(0,nrow=k_regime,ncol=k_regime-1)
    )
  }

  if (overwrite==TRUE){
    use_rvpa_tmb("HMM_SR",overwrite=TRUE)
  }
  if (!file.exists("HMM_SR.dll")) {
    use_rvpa_tmb("HMM_SR")
  }

  obj = TMB::MakeADFun(tmb_data,parameters,DLL="HMM_SR",inner.control=list(maxit=50000,trace=F),silent=TRUE)
  if (length(obj$par)>length(st)) {
    stop("NOT estimable because k > n (k: parameter number, n: sample size")
  }
  opt = nlminb(obj$par,obj$fn,obj$gr,control = list(trace=10,iter.max=10000,eval.max=10000,sing.tol=1e-20))

  Res = list()

  Res$tmb_data = tmb_data
  Res$p0 = parameters
  Res$obj = obj
  Res$opt = opt

  rep = obj$report(opt$par)
  Res$rep = rep

  alpha = rep[["alpha"]]
  a = exp(alpha)
  b = rep[["beta"]]
  sigma = rep[["sigma"]]
  pars=data.frame(regime=1:length(a),a=a,b=b,sd=sigma)
  Res$pars=pars
  regime_prob = (rep[["r_pred"]])
  colnames(regime_prob) = data$year
  Res$regime_prob=round(t(regime_prob),3)
  regime=Rfast::colMaxs(regime_prob)
  names(regime) = data$year
  Res$regime = regime
  Res$trans_prob = rep[["qij"]]

  Res$loglik <- loglik <- -opt$objective
  Res$k <- k <- length(opt$par)
  Res$AIC <- -2*loglik+2*k
  Res$AICc <- Res$AIC+2*k*(k+1)/(length(data$year)-k-1)
  Res$BIC <- -2*loglik+k*log(length(data$year))
  Res$input <- arglist

  # calculate residuals and predicted values
  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)
  SRFV = Vectorize(SRF,vectorize.args = "x")

  pred = sapply(1:nrow(pars), function(i) {
    regime_prob[i,] *log(SRFV(SRdata$SSB,pars$a[i],pars$b[i]))
  }) %>% rowSums %>% exp()

  pred_to_obs = SRdata %>%
    dplyr::rename(Year=year) %>%
    dplyr::mutate(Regime=factor(regime),
                  Pred=pred) %>%
    dplyr::mutate(resid=log(R/pred))

  Res$pred_to_obs <- pred_to_obs
  Res$resid <- pred_to_obs$resid

  ssb.tmp <- seq(from=0,to=max(SRdata$SSB)*max.ssb.pred,length=100)
  pred2 = sapply(1:nrow(pars), function(i) SRFV(ssb.tmp,pars$a[i],pars$b[i])) %>% as.data.frame
  colnames(pred2) <- pars$regime
  pred2 = pred2 %>% mutate(SSB = ssb.tmp) %>%
    pivot_longer(., cols=-SSB,values_to="R") %>%
    mutate(Regime=factor(name)) %>%
    arrange(Regime,SSB) %>%
    dplyr::select(Regime,SSB,R)

  Res$pred <- pred2

  class(Res) <- "hmm_SR"

  return( Res )
}
ichimomo/frasyr documentation built on Dec. 11, 2023, 12:50 a.m.