BRREWABC

Simple estimation using ABC-SMC

library(BRREWABC)

Model definition

compute_dist <- function(x, ss_obs) {
  ss_sim <- c(x[["alpha"]] + x[["beta"]] + rnorm(1, 0, 0.1),
              x[["alpha"]] * x[["beta"]] + rnorm(1, 0, 0.1))
  dist <- sum((ss_sim - ss_obs)^2)
  return(c(dist))
}

model_list <- list("m1" = compute_dist)

Define prior distribution

prior_dist <- list("m1" = list(c("alpha", "unif", 0, 4),
                               c("beta", "unif", 0, 1)))

Create a reference trajectory

sum_stat_obs <- c(2.0, 0.75)

Run abc smc procedure

res <- abcsmc(model_list = model_list,
              prior_dist = prior_dist,
              ss_obs = sum_stat_obs,
              max_number_of_gen = 15,
              nb_acc_prtcl_per_gen = 2000,
              new_threshold_quantile = 0.8,
              experiment_folderpath = "smpl",
              max_concurrent_jobs = 5,
              verbose = FALSE)
#> Warning in scan(file = file, what = what, sep = sep, quote = quote, dec = dec,
#> : EOF within quoted string
#> Warning in scan(file = file, what = what, sep = sep, quote = quote, dec = dec,
#> : EOF within quoted string

Plot results

all_accepted_particles <- res$particles
all_thresholds <- res$thresholds
plot_abcsmc_res(data = all_accepted_particles, prior = prior_dist,
                filename = "smpl/res/figs/smpl_pairplot_all.png", colorpal = "YlGnBu")
#> [1] "Plot saved as 'png'."
plot_thresholds(data = all_thresholds, nb_threshold = 1,
                filename = "smpl/res/figs/smpl_thresholds.png", colorpal = "YlGnBu")
#> [1] "Plot saved as 'png'."
plot_ess(data = all_accepted_particles,filename = "smpl/res/figs/smpl_ess.png", colorpal = "YlGnBu")
#> [1] "Plot saved as 'png'."
#>    gen      ess
#> 1    1 2000.000
#> 2    2 1885.294
#> 3    3 1897.921
#> 4    4 1924.354
#> 5    5 1913.621
#> 6    6 1926.283
#> 7    7 1930.578
#> 8    8 1924.336
#> 9    9 1915.895
#> 10  10 1907.674
#> 11  11 1894.759
#> 12  12 1897.465
#> 13  13 1880.546
#> 14  14 1869.768
#> 15  15 1843.684
plot_densityridges(data = all_accepted_particles, prior = prior_dist,
                   filename = "smpl/res/figs/smpl_densityridges.png", colorpal = "YlGnBu")
#> [1] "Plot saved as 'png'."

Pairplot of all iterations

Threshold evolution over iterations

ESS evolution over iterations

Density estimates for alpha

Density estimates for beta