R/get.parsamp.timedeppar.R

In timedeppar: Infer Constant and Stochastic, Time-Dependent Model Parameters

#' Get a sample of lists of constant and time-dependent parameters from inference results of \code{infer.timedeppar}
#' 
#' This function produces a sample of parameter sets for past and potentially future time points based on the results of
#' class \code{timedeppar} generated by Bayesian inference with the function \code{infer.timedeppar}.
#' For time points used for inference, the sample is a sub-samble of the Markov chain, for future time points of 
#' time-dependent parameters it is a random sample based on the corresponding Ornstein-Uhlenbeck parameters and 
#' constrained at there initial point to the end point of the sub-sample. 
#' 
#' @param  x                  results from the function \code{infer.timedeppar} of class \code{timedeppar}.
#' @param  samp.size          size of the produced sample constructed from the Markov chain stored in the
#'                            object of class \code{timedeppar} omitting the adaptation and burnin phases.
#' @param  n.burnin           number of Markov chain points to omit for density and pairs plots
#'                            (number of omitted points is max(control$n.adapt,n.burnin)).
#' @param  times.new          vector of time points to predict for.
#'                            If no time points are provided, sampling is only from the inference Markov chain;
#'                            if time points are provided, they need to be increasing and start with a larger value 
#'                            than the time points used for inference.
#'                            In the latter case, time-dependent parameters are sampled for the future points and
#'                            appended to the inferred part of the time-dependend parameter.
#' 
#' @return list of\cr
#'         \code{param.maxpost}:   list of constant and time-dependent parameters corresponding to the maximum posterior
#'                                 solution for inference (no extrapolation to the future).\cr
#'         \code{param.maxlikeli}: list of constant and time-dependent parameters corresponding to the solution with
#'                                 maximum observation likelihood found so far.\cr
#'         \code{param.list}:      list of length \code{samp.size} containing lists of constant and time-dependent parameters;
#'                                 for time-dependent parameters sub-sample of the Markov chain for past time points, 
#'                                 sample from Ornstein-Uhlenbeck processes conditioned at the initial point for future
#'                                 time points (see argument \code{times.new}).\cr
#'         \code{param.const}:     sub-sample of constant parameters.\cr
#'         \code{param.timedep}:   list of sub-samples of time-dependent parameters.\cr
#'         \code{param.ou}:        sub-sample of Ornstein-Uhlebeck parameters of the time-dependent parameter(s).\cr
#'         \code{ind.timedeppar}:  indices of time-dependent parameters in the parameter lists.\cr
#'         \code{ind.sample}:      indices of the stored, thinned sample defining the sub-sample.\cr
#'         \code{ind.chain}:       indices of the original non-thinned Markov chain defining the sub-sample.\cr
#'         \code{dot.args}:        ... arguments passed to \code{\link{infer.timedeppar}}; to be re-used for new
#'                                 model evaluations.\cr


get.parsamp <- function(x,samp.size=1000,n.burnin=0,times.new=numeric(0)) UseMethod("get.parsamp")

get.parsamp.timedeppar <- function(x,samp.size=1000,n.burnin=0,times.new=numeric(0))
{
  res <- x
  
  # check input:
  
  if ( ! inherits(res,"timedeppar") )
  {
    warning("get.parsamp.timedeppar: the first argument must be of class timedeppar")
    return(NA)
  }
  errmsg <- character(0)
  if ( is.null(res$param.ini) | 
       is.null(res$sample.param.const) | 
       is.null(res$sample.param.timedep) | 
       is.null(res$sample.param.ou) |
       is.null(res$sample.logpdf) |
       is.null(res$control) )
  {
    errmsg <- c(errmsg,
                 paste("get.parsamp: the first argument must contain objects param.ini, sample.param.const, sample.param.timedep, sample.param.ou, sample.logpdf and control"))
  }
  if ( length(times.new) > 0 )
  {
    if ( any(diff(times.new) <= 0) )
    {
      errmsg <- c(errmsg,
                   "*** get.parsamp: the time points \"times.new\" must be strictly increasing")
    }
  }
  
  # calculate subsample indices:
  
  ind.min <- ceiling(max(res$control$n.adapt,n.burnin)/res$control$thin)
  ind.max <- max(nrow(res$sample.param.const),nrow(res$sample.param.ou))
  if ( ind.max < 1 )
  {
    errmsg <- c(errmsg,
                "*** get.parsamp: no sample data found")
  }
  if ( ind.max < ind.min )
  {
    ind.min <- 1
    warning("get.parsamp.timedeppar: chain elements during adaptation and burnin used")
  }
  if ( ind.max - ind.min + 1 < samp.size )
  {
    samp.size <- ind.max - ind.min + 1
    warning(paste("get.parsamp.timedeppar: sample size reduced to",samp.size))
  }
  if ( length(errmsg) > 0 )
  {
    for ( i in 1:length(errmsg) ) warning(errmsg[i])
    return(NULL)
  }
  ind.sample <- (ind.min:ind.max)[round(seq(from=1,to=ind.max-ind.min+1,length=samp.size))]
  
  # extract maximum posterior and maximum likelihood parameters:
  
  param.all.maxpost <- get.param(res)
  param.all.maxlikeli <- get.param(res,ind.sample=which.max(res$sample.logpdf[,"loglikeliobs"]))
  
  ind.timedeppar <- param.all.maxpost$ind.timedeppar
  
  # extract constant parameter subsample:
  
  param.const <- res$sample.param.const[ind.sample,]
  
  # extract and predict time-dependent parameters to parameter list:
  
  param.list <- list()
  param.ou <- matrix(nrow=samp.size,ncol=3*length(ind.timedeppar))
  for ( i in 1:length(ind.sample) )
  {
    # extract parameters from Markov chain:
    
    param.all.i <- get.param(res,ind.sample=ind.sample[i])
    param.list[[i]]  <- param.all.i$param
    
    # extract OU parameters and predict for new time points:
    
    if ( length(ind.timedeppar) > 0 )
    {
      cn <- character(0)
      for ( j in 1:length(ind.timedeppar) )
      {
        # extract OU parameters:
        
        if ( i == 1 ) cn <- c(cn,paste(names(param.list[[i]])[ind.timedeppar[j]],c("mean","sd","gamma"),sep="_"))
        param.ou[i,(j-1)*3+1] <- param.all.i$param.ou[j,"mean"]
        param.ou[i,(j-1)*3+2] <- param.all.i$param.ou[j,"sd"]
        param.ou[i,(j-1)*3+3] <- param.all.i$param.ou[j,"gamma"]
        
        if ( length(times.new) > 0 )
        {
          # predict for new time points:

          if ( i == 1 )  # check time points
          {
            if ( max(param.list[[1]][[ind.timedeppar[j]]][,1]) > min(times.new) )
            {
              warning("get.parsamp.timedeppar: time points in argument times.new must be larger then time points of inferred parameters")
              return(NA)
            }
          }
          logarithmic <- FALSE
          if ( names(param.list[[i]])[ind.timedeppar[j]] %in% names(res$param.log) )
          {
            if ( res$param.log[names(param.list[[i]])[ind.timedeppar[j]]] ) logarithmic <- TRUE
          }
          val.new <- randOU(mean=param.ou[i,(j-1)*3+1],sd=param.ou[i,(j-1)*3+2],gamma=param.ou[i,(j-1)*3+3],
                            t=c(param.list[[i]][[ind.timedeppar[j]]][nrow(param.list[[i]][[ind.timedeppar[j]]]),1],times.new),
                            yini=param.list[[i]][[ind.timedeppar[j]]][nrow(param.list[[i]][[ind.timedeppar[j]]]),2],log=logarithmic)
          param.list[[i]][[ind.timedeppar[j]]] <- rbind(param.list[[i]][[ind.timedeppar[j]]],as.matrix(val.new[-1,]))
        }
      }
      if ( i == 1 ) colnames(param.ou) <- cn
    }
  }
  
  # copy time-dependent parameters to matrices:
  
  param.timedep <- list()
  if ( length(ind.timedeppar) > 0 )
  {
    for ( j in 1:length(ind.timedeppar) )
    {
      param.timedep[[j]] <- matrix(NA,nrow=1+length(ind.sample),ncol=nrow(param.list[[1]][[ind.timedeppar[j]]]))
      param.timedep[[j]][1,] <- param.list[[1]][[ind.timedeppar[j]]][,1]
      for ( i in 1:length(ind.sample) ) param.timedep[[j]][1+i,] <- param.list[[i]][[ind.timedeppar[j]]][,2]
    }
    names(param.timedep) <- names(res$sample.param.timedep)
  }
  
  dot.args <- list()
  if ( !is.null(res$dot.args) ) dot.args <- res$dot.args
    
  return(list(param.maxpost   = param.all.maxpost$param,
              param.maxlikeli = param.all.maxlikeli$param,
              param.list      = param.list,
              param.const     = param.const,
              param.timedep   = param.timedep,
              param.ou        = param.ou,
              ind.timedeppar  = ind.timedeppar,
              ind.sample      = ind.sample,
              ind.chain       = (ind.sample-1)*res$control$thin,
              dot.args        = dot.args))
}