R/paleoStasis.R

In mcullan/ecic:

#' @export
ecicModel.paleoStasis <- function(model.name, ID = model.name, fix = list()) {
  parameter.names <- c("ns", "theta", "omega", 'vp', 'nn', 'tt') # Define parameters for model.

  #fix = c(fix, list("ns"="ns", "nn"="nn", "vp"="vp","tt"="tt"))
  # Error handling for fixed parameters.
  if(length(fix) > 0){
    # Check if fixed parameters have correct names.
    for (parameter in names(fix) ){
      if ( ! ( parameter %in% parameter.names ) ) {
        message(paste("Model", model.name, "has parameters",
                      paste(parameter.names, collapse = ", "),
                      "so the argument", parameter, "was ignored."))
        fix$parameter <- NULL
      }
    }
    # Assert that vstep is a single numeric value greater than zero.
    if ("vstep" %in% names(fix)){
      assert_that(is.numeric(fix$vstep), length(fix$vstep) == 1, fix$vstep > 0)
    }
    # Assert that mean is a single numeric value.
    if("mstep" %in% names(fix)){
      assert_that(is.numeric(fix$mstep), length(fix$mstep) == 1)
    }
    # Assert that anc is a single numeric value.
    if("anc" %in% names(fix)){
      assert_that(is.numeric(fix$anc), length(fix$anc) == 1)
    }
  } # End error handling for fixed parameters.

  k <- 6 - length(fix)

  return(structure(list(ID = ID, name = "paleoStasis", parameter.names = parameter.names,
                        fixed.parameters = fix,
                        k = k, data.type = "paleoTS"), class = c("ecicModel", "paleoStasis")))
}

#' @export
EstimateParametersMulti.paleoStasis <- function(model, data, ...){
  lapply(data, function(x) EstimateParameters(model, x))
}


#' @export
GenerateData.paleoStasis <- function(n, model, parameters){
  # Computes the log-likelihood and fitted parameters for a dataset and a model.
  #
  # Args:
  #   data:  compatible with the specified model.
  #   n: Sample size.
  #   model: A string specifying the model to compute the log-likehood for.
  #   compress: A boolean specifying if output should be of list or vector type
  #
  # Returns:
  #   A vector/matrix data sample generated by the given model.

  if(n!= parameters$ns){
    stop(paste("n does not match ns parameters for paleoTS object"))
  }
  ns = parameters$ns
  theta = parameters$theta
  omega = parameters$omega
  vp = parameters$vp
  nn = parameters$nn
  tt = parameters$tt
  if (omega < 0 ) omega = 0.0001
  out = sim.Stasis(ns=ns, theta=theta,omega=omega, vp=vp, nn=nn, tt=tt)
  return(out)
}

#' @export
EstimateParameters.paleoStasis = function(model, data){
  if(class(data)!="paleoTS"){
    stop(paste("Data is not of class paleoTS"))
  }
  # if(length(data$mm) != model$nn){
  #   stop(paste("Data is not correct length for ", model$ID,
  #              ", expecting length ", n, ", but input is length ",
  #              length(data), ".", sep = ""))
  # }
  pop.var= pool.var(data)
  fitStasis <- NULL
  attempt <- 1
  while( is.null(fitStasis) && attempt <= 5 ) {
    attempt <- attempt + 1
    tryCatch({
      #seed = sample(1:10000, 1)
      #set.seed(seed)
      fitStasis<- fitSimple(data, "Stasis", pool = FALSE)},
      error = function(e) {
        message("used MLE")
        fitStasis = mle.Stasis(data)}
    )
  }
  paramStasis = fitStasis$parameters
  parameters = list()
  parameters$theta = unname(paramStasis['theta'])
  parameters$omega = unname(paramStasis['omega'])
  parameters$vp = pop.var
  parameters$nn = data$nn
  parameters$ns = length(data$mm)
  parameters$tt = data$tt
  return(parameters)
}

#' @export
logLik.paleoStasis <- function(model, data, compress = F ){
  if(class(data)!="paleoTS"){
    stop(paste("Data is not of class paleoTS"))
  }
  # if(length(data$mm) != model$nn){
  #   stop(paste("Data is not correct length for ", model$ID,
  #              ", expecting length ", n, ", but input is length ",
  #              length(data), ".", sep = ""))
  # }
  pop.var= pool.var(data)
  fitStasis <- NULL
  attempt <- 1
  while( is.null(fitStasis) && attempt <= 5 ) {
    attempt <- attempt + 1
    tryCatch({
      #seed = sample(1:10000, 1)
      #set.seed(seed)
      fitStasis<- fitSimple(data, "Stasis", pool = FALSE)},
      error = function(e) {
        message("used MLE")
        fitStasis = mle.Stasis(data)}
    )
  }
  paramStasis = fitStasis$parameters
  parameters = list()
  parameters$theta = unname(paramStasis['theta'])
  parameters$omega = unname(paramStasis['omega'])
  parameters$vp = pop.var
  parameters$nn = data$nn
  parameters$ns = length(data$mm)
  parameters$tt = data$tt
  logl = fitStasis$logL
  out <- list(log.likelihood = logl, parameters=parameters)
  #out <- c(out, parameters)
  return(out)
}

#' @export
logLikMulti.paleoStasis <- function(model, data, compress = F ){
  if(class(data[[1]])!="paleoTS"){
    stop(paste("Data is not of class paleoTS"))
  }
  # if(length(data$mm) != model$nn){
  #   stop(paste("Data is not correct length for ", model$ID,
  #              ", expecting length ", n, ", but input is length ",
  #              length(data), ".", sep = ""))
  # }
  out = lapply(data, function(x) try(logLik.paleoStasis(model, x)))
  return(out)
}


#' @export
GenerateDataMulti.paleoStasis <- function(n, model, parameters, N){
  lapply(1:N, function(x) GenerateData.paleoStasis(n, model, parameters))
}