R/abc2par.R

In snwikaij/NEMO: What the Package Does (One Line, Title Case)

#' @title Approximate Bayesian algorithm for two parameter models
#'
#' @description Computes the posterior with observations and prior without a likelihood function via an ABC-rejection algorithm.
#'
#' @param observations A numeric vector with c(n, par1, par2) where n denotes the number of observations, par1 the location parameters (i.e. sample mean) and par2 the scale parameter (i.e. sample standard deviation).
#' @param model A specified data generating model with two parameters.
#' @param prior A list with two priors for location an scale parameters.
#' @param nsim Number of simulation by the data generating model.
#' @param qtol The value for the quantile tolerance. A lower value selects simulated parameters that fall closer to the observed parameters. Default is 0.01 and indicates that the closest 1 percent is accepted and the rest rejected.
#' @param adjustment The adjustment method to act if the simulated values actually originated closer from the prior. The current method is the Linear Method (LM), but will also be accompanied by a non-linear RF model.
#' @param print.progress Prints the progress of the number of simulations (nsim) completed. By default is TRUE
#' @param timing Prints the time that was needed to complete the simulations. By default is TRUE.
#' @param seed Default is the Devil. Devils seed.
#'
#' @return A data frame with two four to six vectors.
#'
#' @export abc2par
#'
#' @examples
#' \dontrun{
#'#Create dummy data
#' set.seed(666)
#' x      <- rnorm(5, 50, 5)
#' mxlog  <- mean(log(x))
#' sdxlog <- sd(log(x))
#' mux    <- mean(x)
#' sdx    <- sd(x)
#' n      <- length(x)
#' nsim   <- 200000
#'
#'#Store observed parameters
#' obs <- c(n, mux, sdx)
#'
#Create the data generating model
#' model <- function(n, par1, par2){
#'
#'   simx    <- rnorm(n, par1, par2)
#'
#'   mux     <- mean(simx)
#'
#'   sdx     <- sd(simx)
#'
#'   return(c(mux, sdx))}
#'
#'#Create the priors
#' set.seed(666)
#' priors <- list(prior1  = rlnorm(nsim, 3.4, 0.2),
#'                prior2  = rgamma(nsim, sdx))
#'
#'#Control the priors by eye-balling
#' par(mfrow=c(1,2))
#' hist(priors[[1]], xlab = "", main="Prior1")
#' hist(priors[[2]], xlab = "", main="Prior2")
#' dev.off()
#'
#'#Run the simulations
#' parameters <- abc2par(observations=obs,
#'                       model = model,
#'                       prior = priors,
#'                       nsim = nsim)
#'
#'#Plot posterior mean and sd
#' par(mfrow=c(1,2))
#' hist(parameters$LMadjusted1, main="Posterior mean", xlab="", breaks=15)
#' hist(parameters$LMadjusted2, main="Posterior sd", xlab="", ylab = "", breaks=15)
#'
#'#ABC-result Credible intervals for mean
#' quantile(parameters$LMadjusted1, c(.025, .975))
#'
#'#Confidence intervals for mean
#'c(mux-1.96*sd(x)/sqrt(n), mux+1.96*sd(x)/sqrt(n))}
abc2par <- function(observations, model, prior=NULL, nsim=100000,
                    qtol=0.01, adjustment="LM", print.progress = T,
                    timing = T, seed=666){

  #Start timing
  if(timing == T){stime <- Sys.time()}
  if(is.null(prior)){"Prior cannot be empty"}
  if(is.null(observations)){"Observations cannot be empty"}
  if(length(formalArgs(model)) != length(observations)){"Number of arguments in the model needs to match number of observations"}

  n            <- observations[1]
  lp           <- length(prior)
  simlist      <- list()

  for(i in 1:lp){
    simlist[[i]] <- rep(NA, nsim)}

  if(is.numeric(seed)==T){set.seed(seed)}

  for (p in 1:nsim){

    if(print.progress == T){print(p)}

    sim           <- model(n, prior[[1]][p], prior[[2]][p])

    extsim <- lapply(1:lp, function(x) sim[x])
    for(s in 1:lp){
      simlist[[s]][p] <- extsim[[s]]}

  }

  #Statistics
  statistics <- as.data.frame(matrix(ncol=lp*2, nrow=nsim))
  for(i in 1:lp){
    statistics[,i]       <- simlist[[i]]
    statistics[,i+lp]    <- priors[[i]]}

  colnames(statistics) <- c(paste0(rep("simulate",lp),1:lp), paste0(rep("prior",lp),1:lp))

  #Distance (Epsilon)
  statistics$distance   <- rowSums(sqrt(sweep(statistics[,1:lp], lp, observations[-1])^2))###not yet removed N figured out

  #Select closest quantile qtol
  parameters             <- statistics[order(statistics$distance),][1:c(nsim*qtol),]
  rownames(parameters)   <- NULL

  #Linear adjustment
  if(adjustment == "LM"){
    parameters                 <- cbind.data.frame(as.data.frame(matrix(ncol=lp, nrow=nrow(parameters))), parameters)
    colnames(parameters)[1:lp] <- paste0(rep("LMadjusted",lp),1:lp)
    nr <- 1
    for(c in 1:lp){
      si <- c+lp
      pr <- c+lp+lp
      parameters[,nr]   <- mean(parameters[,pr])+resid(lm(parameters[,pr]~parameters[,si]))
      nr <-nr+1}}

  #Stop timing
  if(timing == T){etime <- Sys.time()
  print(etime-stime)}

  return(parameters)

}