abc2par: Approximate Bayesian algorithm for two parameter models

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

Approximate Bayesian algorithm for two parameter models

Description

Computes the posterior with observations and prior without a likelihood function via an ABC-rejection algorithm.

Usage

abc2par(
  observations,
  model,
  prior = NULL,
  nsim = 100000,
  qtol = 0.01,
  adjustment = "LM",
  print.progress = T,
  timing = T,
  seed = 666
)

Arguments

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).
model	A specified data generating model with two parameters.
prior	A list with two priors for location an scale parameters.
nsim	Number of simulation by the data generating model.
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.
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.
print.progress	Prints the progress of the number of simulations (nsim) completed. By default is TRUE
timing	Prints the time that was needed to complete the simulations. By default is TRUE.
seed	Default is the Devil. Devils seed.

Value

A data frame with two four to six vectors.

Examples

## Not run: 
#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)

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))
## End(Not run)

snwikaij/NEMO documentation built on March 18, 2022, 12:03 a.m.