R/bootstrap.R

In david-borchers/hmltm: Line transect estimation with hidden Markov availability models

#                          ------------------------------
#-------------------------- Start Bootstrap functions ------------------------

#' @title Summarise bootstrap results.
#'
#' @description
#' Produces summarry of bootstrap results
#' 
#' @param bests output from \code{\link{bs.hmltm}}.
#' @param ests hmltm object with points estimates.
#' @param cilevel confidence level for confidence intervals by percentile method. 
#' @param write.csvs if TRUE, writes each output (see Value below) to separate .csv file.
#' @param dir directory to which to write outputs if \code{write.csvs} is TRUE. 
#' @param Dstrat vector containing the names of strata in which density is to be estimated 
#' (must be a subset of \code{unique(dat$stratum)}.)
#'
#' @return If \code{ests} is NULL, returns a list with elements as follows for each statistic 
#' in \code{bests}:
#' \itemize{
#'  \item{nbad} {number of bad estimates of the statistic in question that were excluded from the 
#'  summary. (Bad estimates occur, for example, when a bootstrap sample involves no detections and
#'  the estimation function tries to calculate mean group size.)}
#'  \item{means} {Bootstrap means of the statistic in question.}
#'  \item{cv} {Percentage CVs of the statistic in question.}
#'  \item{se} {SEs of the statistic in question.}
#'  \item{lower} {Lower \code{cilevel} percentiles of the statistic in question.}
#'  \item{upper} {Upper \code{cilevel} percentiles of the statistic in question.}
#' }
#' If \code{ests} is an object of class 'hmltm', returns a data frame with rows only for every 
#' stratum with detections (and the total) and columns as follows:
#' \itemize{
#' \item{Stratum} stratum number
#' \item{n} original number of detections
#' \item{n.L} original encounter rate
#' \item{CV.n.L} percentage coefficient of variation of encounter rate
#' \item{N.grp} original estimate of group abundance
#' \item{CV.N.grp} percentage coefficient of variation of group abundance
#' \item{N.grp.lo} lower bound of group abundance confidence interval
#' \item{N.grp.hi} upper bound of group abundance confidence interval
#' \item{Es} original mean group size estimate
#' \item{CV.Es} percentage coefficient of variation of mean group size
#' \item{Es.lo} lower bound of mean group size confidence interval
#' \item{Es.hi} upper bound of mean group size confidence interval
#' \item{N} original individual abundance estimate
#' \item{CV.N} percentage coefficient of individual abundance
#' \item{N.lo} lower bound of individual abundance confidence interval
#' \item{N.hi} upper bound of individual abundance confidence interval
#'}
#'
#' @export
bootsum <- function(bests,ests=NULL,cilevel=0.95,write.csvs=FALSE,
                    Dstrat=unlist(unique(dimnames(bests)[1])),dir=getwd()){
  if(cilevel<=0 | cilevel>=1) 
    stop("cilevel must be greater than 0 and less than 1.")
  if(!is.null(ests) & !inherits(ests,"hmltm")) 
    stop("ests must be of class 'hmltm'.")
  
  bests = bests[is.element(dimnames(bests)[[1]],Dstrat),,,drop=FALSE] # select only given strata
  
  bdim <- dim(bests)
  B <- bdim[3] # number of bootstraps
  
  # select only strata in Dstrat:
  if(is.element("Total",Dstrat)) {
    warning("Using existing ``Total'' row, not calculating Total from chosen strata.")
    calctotal = FALSE
  } else {
    calctotal = TRUE
    # need to make an array with one more row
    barray = array(rep(NA,(bdim[1]+1)*bdim[2]*bdim[3]),dim=c((bdim[1]+1),bdim[2],bdim[3]))
    for(b in 1:B) barray[,,b] = as.array(rbind(bests[,,b],bests[1,,1]*0))
    dimnames(barray) = list(c(dimnames(bests)[[1]],"Total"),dimnames(bests)[[2]],dimnames(bests)[[3]])
    barray[(bdim[1]+1),,] = apply(barray,c(2,3),sum) # put totals in last rows
    # Density and mean size should not be totals, so correct them here
    barray[(bdim[1]+1),"Dgroups",] = barray[(bdim[1]+1),"Ngroups",]/barray[(bdim[1]+1),"stratum.Area",]
    barray[(bdim[1]+1),"D",] = barray[(bdim[1]+1),"N",]/barray[(bdim[1]+1),"stratum.Area",]
    barray[(bdim[1]+1),"mean.size",] = barray[(bdim[1]+1),"N",]/barray[(bdim[1]+1),"Ngroups",]
    bests = barray
  }
  
  # Keep only bootstrap strata with non-zero sample size:
  ns <- apply(bests[,2,],1,sum) # sum sample sizes by stratum
  ###  keepstrat <- (ns>0) # only consider strata with some detections
  
  if(!is.null(ests)) 
#    keepcols <- c("stratum","n","L","Ngroups","mean.size","N")
    keepcols <- c("n","L","Ngroups","mean.size","N")
  else 
    keepcols <- colnames(bests)
  
  if(!is.null(ests)) 
    bsests <- bests[,keepcols,]
  ###  bsests <- bests[keepstrat,keepcols,]
  else 
    bsests <- bests
  ###  bsests <- bests[keepstrat,,]
  
  # replace L with encounter rate
  Lcol <- which(colnames(bsests[,,1])=="L")
  bsests[,Lcol,] <- bsests[,"n",]/bsests[,"L",] 
  colnames(bsests)[Lcol] <- "n/L"
  
  #  # replace stratum.Area with p
  #  Acol <- which(colnames(bests[,,1])=="stratum.Area")
  #  bsests[,Acol,] <- bsests[,"n",]/(bsests[,"covered.area",]*bsests[,"Dgroups",]) 
  #  colnames(bsests)[Acol] <- "p"
  
  bdim <- dim(bsests)
  nstrat <- bdim[1]
  nests <- bdim[2]
  cv <- matrix(rep(NA,nstrat*nests),ncol=bdim[2])
  rownames <- dimnames(bsests)[[1]]
  rownames[length(rownames)] <- "Total"
  colnames <- dimnames(bsests)[[2]]
  dimnames(cv) <- list(rep("",length(rownames)),colnames)
  nbad <- se <- means <- lower <- upper <- cv
  
  cat("Results from ",B," bootstrap replicates:\n",sep="")
  cat("----------------------------------------\n")
  
  for(i in 1:nstrat) {
    for(j in 1:nests){
      goodests <- na.omit(bsests[i,j,])
      nbad[i,j] <- B-length(goodests)
      means[i,j] <- mean(goodests)
      se[i,j] <- sd(goodests)
      if(means[i,j] == 0) cv[i,j] <- NaN
      else cv[i,j] <- se[i,j]/means[i,j]
      perc <- quantile(goodests,probs = c((1-cilevel)/2,(1-(1-cilevel)/2)))
      lower[i,j] <- perc[1]
      upper[i,j] <- perc[2]
    }
  }
  
  # remove stats of stratum name, and make stratum a character
  nbad <- as.data.frame(nbad,row.names=1:dim(nbad)[1])
  means <- as.data.frame(means,row.names=1:dim(nbad)[1])
  se <- as.data.frame(se,row.names=1:dim(nbad)[1])
  cv <- as.data.frame(cv,row.names=1:dim(nbad)[1])
  lower <- as.data.frame(lower,row.names=1:dim(nbad)[1])
  upper <- as.data.frame(upper,row.names=1:dim(nbad)[1])
  # add stratum column back from rowname:
  nbad = cbind(stratum=rownames,nbad)
  means = cbind(stratum=rownames,means)
  se = cbind(stratum=rownames,se)
  cv = cbind(stratum=rownames,cv)
  lower = cbind(stratum=rownames,lower)
  upper = cbind(stratum=rownames,upper)
#  nbad[,1] <- rownames
#  means[,1] <- rownames
#  se[,1] <- rownames
#  cv[,1] <- rownames
#  lower[,1] <- rownames
#  upper[,1] <- rownames
  
  outsum <- list(nbad=nbad,mean=means,se=se,cv=cv,lower=lower,upper=upper)
  
  if(!is.null(ests)) {
    # Keep only strata that are in bests:
    keepstratnames = row.names(bests)
    keepests = ests$point$ests[is.element(ests$point$ests$stratum,keepstratnames),]
    # need to recalculate totals because ests may have totals of ALL strata:
    totalrow=which(keepests$stratum=="Total")
    if(length(totalrow)==0) { # make total row if it does not exist
      keepests = cbind(keepests,keepests[1,]*0)
      totalrow = dim(keepests)[1]
      keepests$stratum[totalrow] = "Total"
    }
    keepests[totalrow,-1] = apply(keepests[-totalrow,-1],2,"sum") # put totals in last rows; omit 1st col as it is a factor
    # Density and mean size should not be totals, so correct them here
    keepests[totalrow,"Dgroups"] = keepests[totalrow,"Ngroups",]/keepests[totalrow,"stratum.Area"]
    keepests[totalrow,"D"] = keepests[totalrow,"N",]/keepests[totalrow,"stratum.Area"]
    keepests[totalrow,"mean.size"] = keepests[totalrow,"N",]/keepests[totalrow,"Ngroups"]
    outsum <- strat.estable(keepests,outsum$cv)
  }
  
  if(write.csvs){
    dir=as.character(dir)
    wd = getwd()
    if(!file.exists(dir)) dir.create(file.path(wd, dir))
    dir = paste(wd,"/",dir,"/",sep="")
    if(is.null(ests)) {
      write.csv(nbad,file=paste(dir,"nbad.csv",sep=""),row.names=FALSE)
      write.csv(means,file=paste(dir,"means.csv",sep=""),row.names=FALSE)
      write.csv(se,file=paste(dir,"se.csv",sep=""),row.names=FALSE)
      write.csv(cv,file=paste(dir,"cv.csv",sep=""),row.names=FALSE)
      write.csv(lower,file=paste(dir,"lower.csv",sep=""),row.names=FALSE)
      write.csv(upper,file=paste(dir,"upper.csv",sep=""),row.names=FALSE)
    } else write.csv(outsum,file=paste(dir,"bs-summary.csv",sep=""),row.names=FALSE)
  }
  
  return(outsum)
}


#' @title Tabulate bootstrap results.
#'
#' @description
#' Produces brief summarry of bootstrap results
#' 
#' @param est element \code{$point$ests} of output from \code{\link{est.hmltm}}.
#' @param cv \code{$cv} element of output from \code{\link{strat.estable}}. 
#' 
#' @return Returns a data frame with columns as follows:
#' \item{Stratum}{stratum number}
#' \item{n}{original number of detections}
#' \item{n.L}{original encounter rate}
#' \item{CV.n.L}{percentage coefficient of variation of encounter rate}
#' \item{N.grp}{original estimate of group abundance}
#' \item{CV.N.grp}{percentage coefficient of variation of group abundance}
#' \item{N.grp.lo}{lower bound of group abundance confidence interval}
#' \item{N.grp.hi}{upper bound of group abundance confidence interval}
#' \item{E.s}{original mean group size estimate}
#' \item{CV.E.s}{percentage coefficient of variation of mean group size}
#' \item{N}{original individual abundance estimate}
#' \item{CV.N}{percentage coefficient of individual abundance}
#' \item{N.lo}{lower bound of individual abundance confidence interval}
#' \item{N.hi}{upper bound of individual abundance confidence interval}
#' 
#' @export
strat.estable <- function(est,cv){
  dest = dim(est)
  dcv = dim(cv)
  if(dest[1] != dcv[1]) stop(paste("est has ",dest[1]," rows, but cv has ",dcv[1],". They must be the same.",sep=""))
  #  if(dest[2] != dcv[2]) stop(paste("est has ",dest[2]," cols, but cv has ",dcv[2],". They must be the same.",sep=""))
  min <- est[,"n"]
  N <- est[,"N"]
  Ngrp <- est[,"Ngroups"]
  Es <- est[,"mean.size"]
  cv.N <- cv[,"N"]
  cv.Ngp <- cv[,"Ngroups"]
  cv.Es <- cv[,"mean.size"]
  N.ci <- lnci.nmin(min,N,cv.N)
  Ngrp.ci <- lnci.nmin(min,Ngrp,cv.N)
  Es.ci <- lnci.nmin(0,Es,cv.Es)  
  
  outp <- data.frame(Stratum=est$stratum,
                     n=est$n,
                     n.L=signif(est$n/est$L,3),
                     CV.n.L=round(cv[,"n/L"]*100),
                     N.grp=round(est$Ngroups),
                     CV.N.grp=round(cv$Ngroups*100),
                     N.grp.lo=round(Ngrp.ci$lower),
                     N.grp.hi=round(Ngrp.ci$upper),
                     Es=round(est$mean.size,1),
                     CV.Es=round(cv.Es*100),
                     Es.lo=round(Es.ci$lower),
                     Es.hi=round(Es.ci$upper),
                     N=round(est$N),
                     CV.N=round(cv$N*100),
                     N.lo=round(N.ci$lower),
                     N.hi=round(N.ci$upper))
  
  return(outp)
}


#' @title Calculate CV of mean time available.
#'
#' @description
#' Calculates CV of mean time available and unavailable from an hmm.pars object. 
#'
#' @param hmm.pars object of class \code{hmm.pars}. 
#' @param B number of bootstrap replicates to do.
#' 
#' @export
cv.avail <- function(hmm.pars,B=1000){
  n <- dim(hmm.pars$Et)[2]
  b.Et <- array(rep(NA,B*2*n),dim=c(B,2,n),dimnames=list(1:B,State=c("Unavailable","Available"),Animal=1:n))
  cv.a <- cv.u <- rep(NA,n)
  
  for(i in 1:n){
    # get normal parameters corresponding to lognormal(mu,Sigma)
    lN <- Npars.from.lNpars(hmm.pars$Et[,i],hmm.pars$Sigma.Et[,,i])
    
    # resample availability parameters on lognormal scale
    b.Et[,,i] <- exp(mvrnorm(B,lN$mu,lN$Sigma)) 
    
    a <- b.Et[,1,i]/(b.Et[,1,i]+b.Et[,2,i])
    u <- b.Et[,2,i]/(b.Et[,1,i]+b.Et[,2,i])
    cv.a[i] <- cv(a)
    cv.u[i] <- cv(u)
  }
  
  return(list(cv.a=cv.a,cv.u=cv.u))
}

#' @title Bootstrap for hmltm model.
#'
#' @description
#' Stratified nonparameteric bootstrap of hidden Markov line transect model (hmltm) with transects 
#' as the sampling units. 
#'
#' @param hmltm.est output from \code{\link{est.hmltm}}. 
#' @param B number of bootstrap replicates to do.
#' @param hmm.pars.bs output \code{\link{hmmpars.boot}}, containing sets of refitted HMM parameter values.
#' @param bs.trace amount of reporting that \code{\link{optim}} should do while fitting models to
#' bootstrapped datasets. (\code{bs.trace}=0 is no reporting, which is fasters. See component 
#' \code{trace} of parameter \code{control.opt} of function \code{\link{optim}} for more details.)
#' @param report.by frequency with which to report count of number of bootstraps completed.
#' @param fixed.avail whether to treat the hidden Markov availability model as fixed 
#' (\code{fixed.avail}=TRUE) or random (\code{fixed.avail}=FALSE). If \code{fixed.avail}=FALSE, 
#' the availability model is also bootstrapped (see Details below).
#' @param bsfile Name of file to which bootstrap estimates should be written at each iteration 
#' (e.g., in case of crash).
#' 
#' @seealso \code{\link{bootsum}} summarises output from this function.
#' 
#' @details
#' If \code{fixed.avail}=FALSE, then: (1) IF \code{hmltm.est$hmltm.fit$fitpars$hmm.pars$Et} is not 
#' \code{NULL}, the availability process is bootstrapped by drawing pairs of mean times available 
#' (Ea) and unavailable (Eu) from a bivariate lognormal distribution with mean 
#' \code{hmltm.est$hmltm.fit$fitpars$hmm.pars$Et} and standard deviation 
#' \code{hmltm.est$hmltm.fit$fitpars$hmm.pars$Sigma.Et} and converting this to Markov Model transition
#' probability matrix parameters using \code{\link{makePi}}, ELSE IF \code{hmm.pars.bs} is not 
#' \code{NULL}, random samples with replacement, of HMM parameter sets are taken from 
#' \code{hmm.pars.bs}.
#' 
#' @seealso \code{\link{hmmpars.boot}}, which is is used to generate \code{hmm.pars.bs}.
#' 
#' @export
bs.hmltm=function(hmltm.est,B,hmm.pars.bs=NULL,bs.trace=0,report.by=10,fixed.avail=FALSE,bsfile=NULL){
  #######################################################
  ## Produces 3-dim array containing B sets of density ##
  ## and abundance estimates from data extraction      ##
  #######################################################
  
  if(!is.null(bsfile) & !is.character(bsfile)) stop("bsfile must be a character (name of a file)")
  
  dat1 <- hmltm.est$dat # extract original data frame from fitted object
  W.est <- hmltm.est$W.est # extract estimation perp dist truncation
  survey.pars <- hmltm.est$hmltm.fit$fitpars$survey.pars
  hmm.pars <- hmltm.est$hmltm.fit$fitpars$hmm.pars
  control.fit <- hmltm.est$hmltm.fit$fitpars$control.fit
  control.opt <- hmltm.est$hmltm.fit$fitpars$control.opt
  control.opt$trace <- bs.trace
  h.fun <- hmltm.est$hmltm.fit$h.fun
  models <- hmltm.est$hmltm.fit$models
  pars <- hmltm.est$hmltm.fit$fit$par
  
  if(!fixed.avail) {
    # bootstrap availability parameters
    if(!is.null(hmm.pars$Et)){
      cat("Bootstrap with parametric resampling of mean times available and unavailable.\n")
      flush.console()
      
      n <- dim(hmm.pars$Et)[2]
      ns <- bsample(1:n,n,replace=TRUE) # sample animals to use
      b.Et <- array(rep(NA,B*2*n),dim=c(B,2,n),
                    dimnames=list(1:B,State=c("Unavailable","Available"),Animal=1:n))
      
      for(i in ns){
        # get normal parameters corresponding to lognormal(mu,Sigma)
        lN <- Npars.from.lNpars(hmm.pars$Et[,i],hmm.pars$Sigma.Et[,,i]) 
        
        # resample availability parameters on lognormal scale
        b.Et[,,i] <- exp(mvrnorm(B,lN$mu,lN$Sigma)) 
      }
    } else if(!is.null(hmm.pars.bs)) { # resample availability parameters
      cat("Bootstrap with parametric resampling of HMM parameters.\n")
      flush.console()
      
      nhmm <- dim(hmm.pars.bs)[1]
      
      if(nhmm==1) 
        stop(paste("Only one set of hmm pars. need multiple sets of pars if not fixed avail ",
                   "(i.e. if hmm.pars.bs!=NULL)"))
      
      reps <- bsample(1:nhmm,B,replace=TRUE)
    } else {
      warning("No availability parameters to resample from. Treating availability parameters as constant.")
      flush.console()
      fixed.avail <- TRUE
    }
  } else {
    cat("Bootstrap treating HMM parameters as constant.\n")
    flush.console()
  }
  
  # get stuff needed to set up bootstrap datasets
  strat <- unique(dat1$stratum) # unique stratum names
  nrows <- table(dat1$stratum,dat1$transect) # number of rows required for each transect
  get.ntrans <- function(trans) sum(trans>0)
  ntrans <- apply(nrows,1,get.ntrans) # number of transects in each stratum
  nstrat <- length(ntrans) # number of strata
  transects <- as.numeric(colnames(nrows)) # vector of transect numbers
  
  # create bootstrap data, estimate, and store estimates
  nDstrat = length(hmltm.est$point$ests$stratum)-1 # Number of strata for which we want estimates
  Dstrat = hmltm.est$point$ests$stratum[1:nDstrat] # Strata for which we want estimates
  estdim <- dim(hmltm.est$point$ests) #### bitchange
  estdim[2] = estdim[2] - 1 # don't want to keep stratum name column
  bestdim <- c(estdim,B) #### bitchange
#  bestdimnames <- list(as.character(hmltm.est$point$ests$stratum),colnames(hmltm.est$point$ests),1:B)
  bestdimnames <- list(as.character(hmltm.est$point$ests$stratum),colnames(hmltm.est$point$ests)[-1],1:B)
  best <- array(rep(NA,B*prod(estdim)),dim=bestdim,dimnames=bestdimnames)
  bn <- rep(0,B) #### bitchange
  
  for(b in 1:B) { # do bootstraps
    # get hmm pars
    if(!fixed.avail) {
      if(!is.null(hmm.pars$Et)) {
        Pi <- makePi(b.Et[b,1,],b.Et[b,2,]) # make Pi from resampled times
        delta <- matrix(rep(NA,2*n),nrow=2)
        
        for(i in 1:n) 
          delta[,i] <- compdelta(Pi[,,i])
        
        b.hmm.pars <- hmm.pars
        b.hmm.pars$pm <- hmm.pars$pm
        b.hmm.pars$Pi <- Pi
        b.hmm.pars$delta <- delta
        b.hmm.pars$Et <- b.Et[b,,]
      } else {
        b.hmm.pars <- unvectorize.hmmpars(hmm.pars.bs[reps[b],]) # resampled HMM pars
      }
    } else {
      b.hmm.pars <- hmm.pars # fixed availability pars
    }
    
    # to fix naming cock-up when creating hmmpars.bs
    if(is.element("pm",names(b.hmm.pars))) 
      names(b.hmm.pars)[which(names(b.hmm.pars)=="pm")] <- "pm"
    
    # resample transects with replacement
    newtransind <- matrix(rep(NA,nstrat*max(ntrans)),nrow=nstrat) # matrix of indices for resampled transects
    for(st in 1:nstrat) {
      # indices of matrix nrows for resampled transects
      newtransind[st,1:ntrans[st]] <- bsample(which(nrows[st,]>0),ntrans[st],replace=TRUE)
    }
    
    # calc number rows for bootstrap data frame
    newnrows <- 0
    for(st in 1:nstrat) 
      newnrows <- newnrows+sum(na.omit(nrows[st,newtransind[st,]]))
    
    # set up empty bootstrap data frame
    bdat1 <- data.frame(matrix(rep(NA,dim(dat1)[2]*newnrows),nrow=newnrows))
    names(bdat1) <- names(dat1)
    start <- 1
    tno <- 1 # fill in new data frame from old:
    
    for(st in 1:nstrat) 
      for(tr in 1:ntrans[st]) {
        addtrans <- transects[newtransind[st,tr]]
        nadd <- nrows[st,newtransind[st,tr]]
        newi <- start:(start+nadd-1)
        oldi <- which(dat1$stratum==strat[st] & dat1$transect==addtrans)
        bdat1[newi,] <- dat1[oldi,]
        bdat1$transect[newi] <- tno
        start <- start+nadd
        tno <- tno+1
      }
    
    bn[b] <- length(bdat1$object[!is.na(bdat1$object)])
    
    if(b==1) 
      cat("Sample sizes: ")
    
    bhat <- est.hmltm(bdat1,pars=pars,FUN=h.fun,models=models,survey.pars,b.hmm.pars,control.fit,
                      control.opt,notrunc=TRUE,W.est=W.est,Dstrat=Dstrat) 
    if(length(colnames(bhat$point$ests))!=length(colnames(hmltm.est$point$ests)))
      stop("Incompatible columns in bootrapped estimate object and hmltm.est$point$ests; may be using old version of latter?")
    else if(any(colnames(bhat$point$ests)!=colnames(hmltm.est$point$ests)))
      stop("Incompatible columns in bootrapped estimate object and hmltm.est$point$ests; may be using old version of latter?")
    
    for(st in 1:(nDstrat+1)) 
      best[st,,b] <- as.numeric(bhat$point$ests[st,-1]) # exclude stratum column of point$ests
    
    if(b%%report.by==0) {
      cat(paste(bn[b],";",sep=""),"Iterations done:",b,"\n")
      if(b!=B) 
        cat("Sample sizes: ")
    }
    else 
      cat(paste(bn[b],";",sep=""))
  }
  
  class(best) <- c(class(best),"hmltm.bs")
  
  if(!is.null(bsfile)) saveRDS(best,file=bsfile)
  
  return(best)
}


#' @title Normal from logNormal parameters.
#'
#' @description
#' Returns mean and covariance matrix of multivariate normal random variables 
#' which, when logged generate lognormal random variables with mean mu and 
#' covariance matrix Sigma.
#' 
#' @param mu multivariate logNormal distribution mean.
#' @param Sigma multivariate logNormal distribution variace-covarice matrix.
#' 
#' @export
Npars.from.lNpars <- function(mu,Sigma){
  cv2 <- diag(Sigma)/mu^2 # Sigma is variance matrix
  logvar <- log(cv2+1)
  logmu <- log(mu)-logvar/2
  logSigma <- Sigma*0
  np <- length(mu)
  
  for(i in 1:(np-1)) {
    logSigma[i,i] <- logvar[i] # variance of normal
    for(j in (i+1):np) {
      # covariance of normal
      logSigma[i,j] <- log(Sigma[i,j]/(exp(logmu[i]+logmu[j]+(logvar[i]+logvar[j])/2))+1) 
      logSigma[j,i] <- logSigma[i,j]
    } 
  }
  
  logSigma[np,np] <- logvar[np] # variance of normal
  return(list(mu=logmu,Sigma=logSigma))
}


#' @title Parameric resample availability HMM.
#'
#' @description
#' Resamples a hidden Markov model (HMM) data from fitted variance-covariance
#' matrices (one per animal).
#' 
#' @param availhmm list with availability HMM paramters, as for the \code{hmm.pars} parameter of 
#' \code{\link{est.hmltm}}.
#' @param animals a vector of integers indicating which of the elements of \code{adat} are to be
#' resampled (e.g. animals=c(1,3,4) means the time series of animals 1, 3 and 4 only will be used).
#' @param Pi.only If TRUE, resamples transition probability matrix (TPM) 
#' parameters but not state-dependent parameters (which are taken as fixed). The 
#' reason for this is that with availability datasets in which the state-dependent 
#' parameters are close to 1 and 0, the full variance-covariance matrix is often 
#' singular, although the TPM component is not.
#' @param seed random number seed (integer).
#' @param nperow number of iterations printed on one row if printprog=TRUE.
#' @param printprog if TRUE, prints progress through animals as it resamples.
#'
#' @details
#' Resamples HMM paramters assuming asymptotic normality of parameters. 
#' Constructs a new hmm.pars object from the fitted HMM parameters.
#' 
#' @return A list with the same format as the input \code{availhmm}
#' 
#' @export
parametric.resample.hmmpars <- function(availhmm,animals,Pi.only=TRUE,seed=NULL,nperow=10,printprog=TRUE){
  b.availhmm <- availhmm[1:3] # initialise bootstrap HMM parameters (excluding vcv component)
  
  # filter out animals not chosen
  nstates = dim(b.availhmm$Pi[,,1])[1]
  b.availhmm$Pi <- b.availhmm$Pi[,,animals]
  b.availhmm$pm <- b.availhmm$pm[,animals]
  b.availhmm$delta <- b.availhmm$delta[,animals]
  
  # simulate and fit new HMMs
  newanimal <- 1
  for(animal in animals) {
    
    if(Pi.only) {
      mu = c(availhmm$Pi[1,1,animal],availhmm$Pi[2,1,animal])
      vcv = availhmm$vcv[1:nstates,1:nstates,animal]
      b.pars = mvrnorm(n=1,mu=logit(mu),Sigma=vcv)
      b.availhmm$pm[,newanimal] <- availhmm$pm[,animal]
      b.availhmm$Pi[,,newanimal] <- matrix(c(inv.logit(b.pars),1-inv.logit(b.pars)),nrow=nstates)
      b.availhmm$delta[,newanimal] <- compdelta(b.availhmm$Pi[,,newanimal])
    } else {
      mu = c(availhmm$Pi[1,1,animal],availhmm$Pi[2,1,animal],availhmm$pm[,animal])
      vcv = availhmm$vcv[,,animal]
      b.pars = mvrnorm(n=1,mu=logit(mu),Sigma=vcv)
      b.availhmm$pm[,newanimal] <- inv.logit(b.pars[(nstates+1):length(b.pars)])
      b.availhmm$Pi[,,newanimal] <- matrix(c(inv.logit(b.pars[1:nstates]),1-inv.logit(b.pars[1:nstates])),nrow=nstates)
      b.availhmm$delta[,newanimal] <- compdelta(b.availhmm$Pi[,,newanimal])
    }
    if(printprog) {
      if(animal==1) 
        cat("Individuals done:")
      if((animal%/%nperow)*nperow==animal) # got an exact multiple of 10
        if(animal>1) 
          cat(": Total=",newanimal,"\nIndividuals done:")
    }
    
    cat(" *")
    newanimal <- newanimal+1
  }
  
  cat("\n Total=",newanimal," individuals done\n")
  return(b.availhmm)  
}


#' @title Resample availability HMM.
#'
#' @description
#' Resamples a hidden Markov model (HMM) data from multiple observed availability time series 
#' (one per animal).
#' 
#' @param availhmm list with availability HMM paramters, as for the \code{hmm.pars} parameter of 
#' \code{\link{est.hmltm}}.
#' @param adat list of availability data time series. The ith element of the list must be named 
#' \code{$ai} and must be a vector of 0s and 1s, comprising the ith animal's time series, with 0 
#' corresponding to being unavailable and 1 to being available.
#' @param animals a vector of integers indicating which of the elements of \code{adat} are to be
#' resampled (e.g. animals=c(1,3,4) means the time series of animals 1, 3 and 4 only will be used).
#' @param seed random number seed (integer).
#' @param nperow number of iterations printed on one row if printprog=TRUE.
#' @param printprog if TRUE, prints progress through animals as it resamples.
#'
#' @details
#' Simulates a new series of availability observations (0s and 1s) using functions \code{dthmm} 
#' and \code{simulate} from library \code{\link{HiddenMarkov}}, then fits a HMM to these data
#' using function \code{BaumWelch} from library \code{\link{HiddenMarkov}}. Constructs a new 
#' hmm.pars object from the fitted HMM parameters.
#' 
#' @return A list with the same format as the input \code{availhmm}
#' 
#' @examples
#' data(bowhead.hmm.pars)
#' data(bowhead.adat)
#' animals=1:8
#' resamp=resample.hmmpars(bowhead.hmm.pars,bowhead.adat,animals)
#' 
#' @importFrom HiddenMarkov dthmm
#' @importFrom HiddenMarkov BaumWelch
#' 
#' @export
resample.hmmpars <- function(availhmm,adat,animals,seed=NULL,nperow=10,printprog=TRUE){
  b.availhmm <- availhmm # initialise bootstrap HMM parameters
  
  # filter out animals not chosen
  b.availhmm$Pi <- b.availhmm$Pi[,,animals]
  b.availhmm$pm <- b.availhmm$pm[,animals]
  b.availhmm$delta <- b.availhmm$delta[,animals]
  
  # simulate and fit new HMMs
  newanimal <- 1
  for(animal in animals) {
    hmmobj <- dthmm(adat[[animal]],availhmm$Pi[,,animal],availhmm$delta[,animal],"binom",
                    pm=list(prob=availhmm$pm[,animal]),pn=list(size=rep(1,length(adat[[animal]]))),
                    discrete = TRUE)
    
    sdat <- simulate(hmmobj, nsim=length(hmmobj$x),seed=seed)
    fitanimal <- BaumWelch(sdat,control=list(maxiter=500,tol=1e-05,prt=FALSE,posdiff=TRUE,
                                             converge = expression(diff < tol)))
    
    b.availhmm$pm[,newanimal] <- fitanimal$pm$prob
    b.availhmm$Pi[,,newanimal] <- fitanimal$Pi
    b.availhmm$delta[,newanimal] <- compdelta(fitanimal$Pi)
    
    if(printprog) {
      if(animal==1) 
        cat("Individuals done:")
      if((animal%/%nperow)*nperow==animal) # got an exact multiple of 10
        if(animal>1) 
          cat(": Total=",animal,"\nIndividuals done:")
    }
    
    cat(" *")
    newanimal <- newanimal+1
  }
  
  cat("\n Total=",animal," individuals done\n")
  return(b.availhmm)  
}




#' @title Reformat hmm.pars object as a vector.
#'
#' @description
#' Reformats hmm.pars object as a vector - so that it can easily be written to file.
#' 
#' @param hmmpars list with availability HMM paramters, as for the \code{hmm.pars} parameter of 
#' \code{\link{est.hmltm}}.
#' 
#' @return A numeric vector.
#' 
#' @export
vectorize.hmmpars <- function(hmmpars) {
  Pi <- hmmpars$Pi
  pm <- hmmpars$pm
  delta <- hmmpars$delta
  dims <- dim(Pi)
  
  if(length(dims)!=3) 
    stop("Pi must be a 3-D array")
  if(dims[1]!=dims[2]) 
    stop("First two dimensions of Pi unequal")
  if(dim(pm)[1]!=dims[1] | dim(pm)[2]!=dims[3]) 
    stop("dim(pm) inconsistent with dim(Pi)")
  if(dim(delta)[1]!=dims[1] | dim(delta)[2]!=dims[3]) 
    stop("dim(delta) inconsistent with dim(Pi)")
  
  return(c(as.vector(dim(Pi)),as.vector(Pi),as.vector(pm),as.vector(delta)))
}

#' @title Reformat vector as hmm.pars object.
#'
#' @description
#' Reformats vector that was vectorised using \code{\link{vectorize.hmmpars}}, as a hmm.pars object.
#' 
#' @param hv vector that was vectorised using \code{\link{vectorize.hmmpars}}.
#' 
#' @return A hmm.pars object (a list).
#' 
#' @export
unvectorize.hmmpars <- function(hv) {
  m3d <- hv[1:3]
  m2d <- c(m3d[1],m3d[3])
  m3size <- prod(m3d)
  m2size <- prod(m2d)
  Pi <- array(hv[(3+1):(3+m3size)],dim=m3d)
  pm <- array(hv[(3+m3size+1):(3+m3size+m2size)],dim=m2d)
  delta <- array(hv[(3+m3size+m2size+1):(3+m3size+m2size+m2size)],dim=m2d)
  return(list(pm=pm,Pi=Pi,delta=delta))
}

#' @title Bootstrap availability HMM.
#'
#' @description
#' Bootstraps hidden Markov model (HMM) data from multiple observed availability time series 
#' (one per animal).
#' 
#' @param availhmm list with availability HMM paramters, as for the \code{hmm.pars} parameter of 
#' \code{\link{est.hmltm}}.
#' @param adat list of availability data time series. The ith element of the list must be named 
#' \code{$ai} and must be a vector of 0s and 1s, with 0 corresponding to being unavailable and 1 to 
#' being available.
#' @param animals a vector of integers indicating which of the elements of \code{adat} are to be
#' resampled.
#' @param B number of bootstraps to perform (integer).
#' @param seed random number seed (integer).
#' @param printprog if TRUE prints progress through animals as it resamples.
#'
#' @details
#' Simulates a new series of availability observations (0s and 1s) using functions \code{dthmm} 
#' and \code{simulate} from library \code{\link{HiddenMarkov}}, then fits a HMM to these data
#' using function \code{BaumWelch} from library \code{\link{HiddenMarkov}}. Constructs a new 
#' hmm.pars object from the fitted HMM parameters. 
#' 
#' Does the above \code{B} times, each time reformtting the hmm.pars object as a vector using 
#' \code{\link{vectorize.hmmpars}} and then entering this as the next row in a matrix of dimension 
#' \code{B}xT, where T=\code{3+length(animals)*(nstate^2+nstate*2)} and \code{nstate} is the 
#' number of states in the HMM.
#' 
#' @return A matrix of dimension \code{B}xT, where T is as described above.
#' 
#' @examples
#' data(bowhead.hmm.pars)
#' data(bowhead.adat)
#' animals=1:8
#' bs=hmmpars.boot(bowhead.hmm.pars,bowhead.adat,animals,B=3)
#' 
#' @export
hmmpars.boot <- function(availhmm,adat,animals,seed=NULL,B,printprog=TRUE){
  # initialise matrix of correct dimensions:
  na <- length(adat)
  #  ncol=1
  #  for(i in 1:na) if(length(adat[[i]])>ncol) ncol=length(adat[[i]])
  
  onelength <- length(vectorize.hmmpars(list(Pi=availhmm$Pi[,,1,drop=FALSE],
                                             pm=availhmm$pm[,1,drop=FALSE],
                                             delta=availhmm$delta[,1,drop=FALSE])))
  
  ncol <- 3+(onelength-3)*length(animals) # 3 numbers specify vectorized length parameters
  hmmat <- matrix(rep(NA,ncol*B),ncol=ncol)
  
  # do the bootstrapping
  for(b in 1:B) {
    hmp <- resample.hmmpars(availhmm,adat,animals,seed=NULL)
    hmmvec <- vectorize.hmmpars(hmp)
    #    nobs=length(hmmvec)
    #    hmmat[b,nobs]=hmmvec
    hmmat[b,] <- hmmvec
    
    if(printprog) 
      cat("Resamples done: ",b,"\n")
  }
  
  cat("\n Total=",b," Resamples done\n")
  return(hmmat)
}



#' @title Parameteric bootstrap availability HMM.
#'
#' @description
#' Bootstraps hidden Markov model (HMM) data from estimated parameters and their
#' variance-covariance matrix (one of each for each animal).
#' 
#' @param availhmm list with availability HMM parameters, as for the \code{hmm.pars} parameter of 
#' \code{\link{est.hmltm}}. The list must contain these elements:
#' \itemize{
#' \item{Pi:}{ a 3D matrix of Markov chain transition probability matrices, with 
#' element [,,i] containing the 2D matrix (Pi) for animal i 
#' (state 1=unavailable, state 2=available).}
#' \item{pm:}{ a matrix with element [,i] containing the state-dependent Bernoulli 
#' ``success'' parameters (p) for each state (first state=unavailable, 2nd=available)
#' for animal i.}
#' \item{delta:}{ a matrix with element [,i] containing the stationary disribution
#' of the Markov chain (first state=unavailable, 2nd=available).}
#' \item{vcv:}{ a 3D array with element [,,i] containing the estimated variance-
#' covariance matrix of the following (in this order): logit(Pi[1,1]),
#' logit(Pi[2,1]),logit(p[1]),logit(p[2]).}
#' }
#' @param animals a vector of integers indicating which of the elements of \code{adat} are to be
#' resampled.
#' @param B number of bootstraps to perform (integer).
#' @param seed random number seed (integer).
#' @param printprog if TRUE prints progress through animals as it resamples.
#'
#' @details
#' Resamples HMM paramters assuming asymptotic normality of parameters. Constructs a new 
#' hmm.pars object from the fitted HMM parameters. 
#' 
#' Does the above \code{B} times, each time reformtting the hmm.pars object as a vector using 
#' \code{\link{vectorize.hmmpars}} and then entering this as the next row in a matrix of dimension 
#' \code{B}xT, where T=\code{3+length(animals)*(nstate^2+nstate*2)} and \code{nstate} is the 
#' number of states in the HMM.
#' 
#' @return A matrix of dimension \code{B}xT, where T is as described above.
#' 
#' @export
hmmpars.paramtetric.boot <- function(availhmm,animals,seed=NULL,B,printprog=TRUE){
  # initialise matrix of correct dimensions:
  na <- dim(availhmm$Pi)[3]
  #  ncol=1
  #  for(i in 1:na) if(length(adat[[i]])>ncol) ncol=length(adat[[i]])
  
  onelength <- length(vectorize.hmmpars(list(Pi=availhmm$Pi[,,1,drop=FALSE],
                                             pm=availhmm$pm[,1,drop=FALSE],
                                             delta=availhmm$delta[,1,drop=FALSE])))
  
  ncol <- 3+(onelength-3)*length(animals) # 3 numbers specify vectorized length parameters
  hmmat <- matrix(rep(NA,ncol*B),ncol=ncol)
  
  # do the bootstrapping
  for(b in 1:B) {
    hmp <- parametric.resample.hmmpars(availhmm,animals,seed=NULL)
    hmmvec <- vectorize.hmmpars(hmp)
    #    nobs=length(hmmvec)
    #    hmmat[b,nobs]=hmmvec
    hmmat[b,] <- hmmvec
    
    if(printprog) 
      cat("Resamples done: ",b,"\n")
  }
  
  cat("\n Total=",b," Resamples done\n")
  return(hmmat)
}



#' @title Detection probability bootstrap with availability process times.
#'
#' @description
#' Nonparametric bootstrap of detection data with estimation of detection probabilities. If 
#' \code{fixed.avail}=FALSE, does parametric resampling of mean times available and unavailable for 
#' every resample of detection data, else treats these mean times as fixed.
#' 
#' @param dat detection data frame constructed by removing all rows with no detections from a 
#' data frame of the sort passed to \code{\link{est.hmltm}}.
#' @param pars starting parameter values, as for \code{\link{est.hmltm}}.
#' @param hfun detection hazard function name; same as argument \code{FUN} of \code{\link{est.hmltm}}.
#' @param models detection hazard covariate models, as for \code{\link{est.hmltm}}.
#' @param survey.pars survey parameters, as for \code{\link{est.hmltm}}.
#' @param hmm.pars availability hmm parameters, as for \code{\link{est.hmltm}}. Must have elements
#' \code{$Et} and \code{$Sigma.Et}
#' @param control.fit list controlling fit, as for \code{\link{est.hmltm}}.
#' @param control.opt list controlling function \code{\link{optim}}, as for \code{\link{est.hmltm}}.
#' @param fixed.avail if TRUE, hmm.pars is treated as fixed, else element \code{$Et} is parametrically 
#' resampled.
#' @param B number of bootstrap replicates.
#' 
#' @details
#' The rows of data frame \code{dat} are resampled with replacement to create new data frames with as 
#' many detections as were in \code{dat}. If \code{fixed.avail}=TRUE, then a pair of new mean times 
#' available and unavailable (\code{$Et}s) are generated for each resampled data frame, by resampling 
#' parametrically from a logNormal distribution with mean \code{hmm.pars$Et} and variance-covariance
#' matrix \code{hmm.pars$Sigma.Et}.
#' 
#' Function \code{\link{fit.hmltm}} is called to estimate detection probabilities and related things
#' for every bootstrap resample.
#' 
#' @return 
#' A list with the following elements:
#' \itemize{
#' \item{callist:}{ input reflection: everything passed to the function, bundled into a list}
#' \item{bs:}{ a list containing (a) a Bxn matrix \code{$phats} in which each row is the estimated 
#' detection probabilities for each of the n bootstrapped detections, (b) a Bxn matrix \code{$pars} 
#' in which each row is the estimated detection hazard parameters, (c) the following vectors 
#' of length B with estimates from each bootstrap: \code{$p0} (mean estimated p(0) over all 
#' detections), \code{$phat} (mean estimated detection probability over all detections), and (d) 
#' a Bx2 matrix \code{$b.Et} in which each row is the mean times unavailable and available.
#' }
#' }
#' 
#' @importFrom MASS mvrnorm
#' @importFrom HiddenMarkov compdelta
#' 
#' @export
bootstrap.p.with.Et <- function(dat,pars,hfun,models,survey.pars,hmm.pars,
                                control.fit,control.opt,fixed.avail=FALSE,B=999){
  n <- length(dat$x)
  npar <- length(pars)
  b.p0 <- b.phat <- rep(NA,B)
  b.pars <- matrix(rep(NA,B*npar),ncol=npar)
  b.phats <- matrix(rep(NA,B*n),ncol=n)# matrix for detection probs of each individual.
  
  # bootstrap availability parameters
  if(!fixed.avail) {
    # get normal parameters corresponding to lognormal(mu,Sigma)
    lN <- Npars.from.lNpars(hmm.pars$Et,hmm.pars$Sigma.Et) 
    
    # resample availability parameters on lognormal scale
    b.Et <- exp(mvrnorm(B,lN$mu,lN$Sigma)) 
  }
  
  # resample sightings data and re-estimate, using a resample of availability paramters:
  for(nb in 1:B){
    # resample detection locations
    samp.ind <- bsample(1:n,size=n,replace=TRUE) # resample sightings data indices with replacement
    b.dat <- dat[samp.ind,]# get resampled data
    
    # create new hmm.pars object with resampled availability parameters
    if(!fixed.avail) {
      pi21 <- 1/b.Et[nb,2]
      pi12 <- 1/b.Et[nb,1]
      Pi <- matrix(c((1-pi12),pi12,pi21,(1-pi21)),nrow=2,byrow=TRUE)
      delta <- compdelta(Pi)
      pm <- c(0.0,1.0)
      b.hmm.pars <- list(pm=pm,Pi=Pi,delta=delta,b.Et[nb],hmm.pars$Sigma.Et)
    }else {
      b.hmm.pars <- hmm.pars
    }
    
    # refit model
    b.fit <- fit.hmltm(b.dat,pars=pars,FUN=hfun,models=models,survey.pars=survey.pars,
                       hmm.pars=b.hmm.pars,control.fit=control.fit,control.optim=control.opt)
    
    b.p0[nb] <- b.fit$pzero
    b.phat[nb] <- b.fit$phat
    b.phats[nb,] <- b.fit$phats
    b.pars[nb,] <- b.fit$fit$par
    
    cat(paste("done",nb,"\n"))
    flush.console()
  }
  
  # package results and return
  callist <- list(dat=dat,pars=pars,hmm.pars=hmm.pars,hfun=hfun,models=models,
                  survey.pars=survey.pars,control.fit=control.fit,control.optim=control.opt)
  bs <- list(phats=b.phats,pars=b.pars,p0=b.p0,phat=b.phat,b.Et=b.Et)
  
  return(list(callist=callist,bs=bs))
}




#' @title Detection probability bootstrap with availability HMM.
#'
#' @description
#' Nonparametric bootstrap of detection data with re-estimation of detection probabilities. If 
#' \code{fixed.avail}=FALSE, does nonparametric resampling of availability HMM parameters contained
#' in hmm.pars.bs for every resample of detection data.
#' 
#' @param dat detection data frame constructed by removing all rows with no detections from a 
#' data frame of the sort passed to \code{\link{est.hmltm}}.
#' @param pars starting parameter values, as for \code{\link{est.hmltm}}.
#' @param hfun detection hazard function name; same as argument \code{FUN} of \code{\link{est.hmltm}}.
#' @param models detection hazard covariate models, as for \code{\link{est.hmltm}}.
#' @param survey.pars survey parameters, as for \code{\link{est.hmltm}}.
#' @param hmm.pars.bs multiple sets of availability hmm parameters, as output by 
#' \code{\link{hmmpars.boot}}, OR a single set of hmm parameters in appropriate format
#' @param control.fit list controlling fit, as for \code{\link{est.hmltm}}.
#' @param control.opt list controlling function \code{\link{optim}}, as for \code{\link{est.hmltm}}.
#' @param fixed.avail if TRUE, hmm.pars is treated as fixed, else element \code{$Et} is parametrically 
#' resampled.
#' @param B number of bootstrap replicates.
#' @param silent argument of function \code{\link{try}}, controlling error message reporting.
#' 
#' @details
#' The rows of data frame \code{dat} are resampled with replacement to create new data frames with as 
#' many detections as were in \code{dat}. If \code{fixed.avail}=TRUE, then a new set of availability 
#' HMM parameters is obtainded by sampling iwth replacement from \code{hmm.pars.bs}.
#' 
#' Function \code{\link{est.hmltm}} is called to estimate detection probabilities and related things
#' for every bootstrap resample.
#' 
#' @return 
#' A list with the following elements:
#' \itemize{
#' \item{callist:}{ input reflection: everything passed to the function, bundled into a list}
#' \item{bs:}{ a list containing (a) a Bxn matrix \code{$phats} in which each row is the estimated 
#' detection probabilities for each of the n bootstrapped detections, (b) a Bxn matrix \code{$pars} 
#' in which each row is the estimated detection hazard parameters, and (c) the following vectors 
#' of length B with estimates from each bootstrap: \code{$Et} (mean times available and unavailable), 
#' \code{$p0} (mean estimated p(0) over all detections), \code{$phat} (mean estimated detection
#' probability over all detections), \code{$convergence} convergence diagnostic from \code{optim}.}
#' }
#' 
#' @export
bootstrap.p.with.hmm <- function(dat,pars,hfun,models,survey.pars,hmm.pars.bs,
                                 control.fit,control.opt,fixed.avail=FALSE,B=999,silent=FALSE){
  n <- length(dat$x)
  npar <- length(pars)
  conv <- b.p0 <- b.phat <- rep(NA,B)
  b.pars <- matrix(rep(NA,B*npar),ncol=npar)
  b.phats <- matrix(rep(NA,B*n),ncol=n) # matrix for detection probs of each individual.
  
  # bootstrap availability parameters
  if(!fixed.avail) { # resample availability parameters
    nhmm <- dim(hmm.pars.bs)[1]
    if(nhmm==1) stop("Only one set of hmm pars. need multiple sets of pars if not fixed avail.")
    reps <- bsample(1:nhmm,B,replace=TRUE)
  }
  for(nb in 1:B) {
    if(!fixed.avail)
      b.hmm.pars <- unvectorize.hmmpars(hmm.pars.bs[reps[nb],])
    else {
      if(is.list(hmm.pars.bs)) b.hmm.pars <- hmm.pars.bs # here if not bootstrapped but just one set of pars
      else b.hmm.pars <- unvectorize.hmmpars(hmm.pars.bs[1,]) # here if bootstrapped sets of pars
    }
    
    # to fix naming cock-up when creating hmmpars.bs
    if(is.element("pm",names(b.hmm.pars))) 
      names(b.hmm.pars)[which(names(b.hmm.pars)=="pm")] <- "pm" 
    
    # resample detection locations
    samp.ind <- bsample(1:n,size=n,replace=TRUE) # resample sightings data indices with replacement
    b.dat <- dat[samp.ind,,drop=FALSE]# get resampled data
    names(b.dat) <- names(dat)
    
    # refit model
    b.fit <- try(fit.hmltm(b.dat,pars=pars,FUN=hfun,models=models,survey.pars=survey.pars,
                           hmm.pars=b.hmm.pars,control.fit=control.fit,control.optim=control.opt),
                 silent=silent)
    
    if((class(b.fit)=="try-error")) {
      conv[nb] <- -999
      b.p0[nb] <- -999
      b.phat[nb] <- -999
      b.phats[nb,] <- rep(-999,n)
      b.pars[nb,] <- rep(-999,length(pars))
    } else {
      conv[nb] <- b.fit$fit$convergence
      b.p0[nb] <- b.fit$p[1]
      b.phat[nb] <- b.fit$phat
      b.phats[nb,] <- b.fit$phats
      b.pars[nb,] <- b.fit$fit$par
      cat(paste("done",nb,"\n"))
    }
    flush.console()
  }
  
  # package results and return
  callist <- list(dat=dat,pars=pars,hmm.pars=hmm.pars.bs,hfun=hfun,survey.pars=survey.pars,
                  control.fit=control.fit,control.optim=control.opt)
  
  bs <- list(phats=b.phats,pars=b.pars,p0=b.p0,phat=b.phat,convergence=conv)
  return(list(callist=callist,bs=bs))
}


#' @title Summarise detection probability bootstrap results.
#'
#' @description
#' Uses bootstrap results from \code{\link{bootstrap.p.with.Et}} or 
#' \code{\link{bootstrap.p.with.hmm}} to work out bootstrap means, variance estimates, CVs and
#' confidence intervals.
#' 
#' @param bs output from \code{\link{bootstrap.p.with.Et}} or \code{\link{bootstrap.p.with.hmm}}.
#' @param probs lower and upper percentile points for confidence interval reporting.
#' @param pcut minimum estimated detection probability to use. This is a quick and dirty method to 
#' robustify against small detection probability estimates skewing the distribution of \code{1/phat} 
#' badly for small samples. It is ad-hoc. If you use it, do histogram of $bs$phat to see if there is 
#' a reasonable cutpoint.
#' 
#' @return 
#' Returns a list with elements
#' \itemize{
#' \item{nboot:}{ number of bootstrap estimates used in constructing bootstrap statistics.}
#' \item{nbad:}{ number of bad estimates excluded from results.}
#' \item{parcov:}{ parameter estimate variance-covariance matrix.}
#' \item{parcorr:}{ parameter estimate correlation matrix.}
#' \item{bests:}{ bootstrap estimate statistics, comprising meand, standard error, percentage CV and 
#' confidence interval limits for: (a) estimated mean detection probability, (a) 1/(estimated mean 
#' detection probability), (c) estimated p(0), ad (d) detection hazard function parameters.}
#' }
#' 
#' @export
bootsum.p <- function(bs,probs=c(0.025,0.975),pcut=0){
  ################################################################################
  ## bs is output from bootstrap.p.with.Et() or bootstrap.with.hmm()            ##
  ## pcut is a quick and dirty min phat to allow - robustifies 1/phat for small ##
  ## samples, although it is ad-hoc. Do hist of $bs$phat to see if there is a   ##
  ## reasonable cutpoint.                                                       ##
  ################################################################################
  
  nboot <- length(bs$bs$p0)
  
  if(is.null(bs$bs$convergence)) 
    keep <- which(bs$bs$p0>=0 & bs$bs$phat>pcut)
  else 
    keep <- which(bs$bs$p0>=0 & bs$bs$convergence==0 & bs$bs$phat>pcut)
  
  nbad <- nboot-length(keep)
  npar <- dim(bs$bs$par)[2]
  cinames <- paste(as.character(probs*100),"%",sep="")  
  bests <- matrix(rep(NA,(3+npar)*5),ncol=5)
  colnames(bests) <- c("mean","std.err.","%CV",cinames)
  parnames <- paste("par",as.character(1:npar),sep="")
  rownames(bests) <- c("1/phat","phat","p(0)",parnames)
  
  # relative density
  bests[1,1] <- mean(1/bs$bs$phat[keep])
  bests[1,2] <- sd(1/bs$bs$phat[keep])
  bests[1,3] <- sd(1/bs$bs$phat[keep])/mean(1/bs$bs$phat[keep])*100
  bests[1,4:5] <- quantile(1/bs$bs$phat[keep],probs=probs)
  
  # detection probability
  bests[2,1] <- mean(bs$bs$phat[keep])
  bests[2,2] <- sd(bs$bs$phat[keep])
  bests[2,3] <- sd(bs$bs$phat[keep])/mean(bs$bs$phat[keep])*100
  bests[2,4:5] <- quantile(bs$bs$phat[keep],probs=probs)
  
  # p(0)
  bests[3,1] <- mean(bs$bs$p0[keep])
  bests[3,2] <- sd(bs$bs$p0[keep])
  bests[3,3] <- sd(bs$bs$p0[keep])/mean(bs$bs$p0[keep])*100
  bests[3,4:5] <- quantile(bs$bs$p0[keep],probs=probs)
  
  # parameters
  for(i in 1:npar) {
    bests[3+i,1] <- sd(bs$bs$par[keep,i])
    bests[3+i,2] <- mean(bs$bs$par[keep,i])
    bests[3+i,3] <- sd(bs$bs$par[keep,i])/mean(bs$bs$par[keep,i])*100
    bests[3+i,4:5] <- quantile(bs$bs$par[keep,i],probs=probs)
  }
  
  parcov <- cov(bs$bs$par)
  parcorr <- cov2cor(parcov)
  rownames(parcov) <- rownames(parcorr) <- colnames(parcov) <- colnames(parcorr) <- parnames
  
  return(list(nboot=nboot,nbad=nbad,bests=bests,parcov=parcov,parcorr=parcorr))
}

#-------------------------- End Bootstrap functions ------------------------
#                          ------------------------------