Using cmdstanr in SimDesign

library(SimDesign)
library(cmdstanr)

[Update: Use parallel computing with two cores.]

Adapted from https://cran.r-project.org/web/packages/SimDesign/vignettes/SimDesign-intro.html

See https://mc-stan.org/cmdstanr/articles/cmdstanr.html for using cmdstanr

Design <- createDesign(sample_size = c(30, 60, 120, 240), 
                       distribution = c('norm', 'chi'))
Design

## # A tibble: 8 × 2
##   sample_size distribution
##         <dbl> <chr>       
## 1          30 norm        
## 2          60 norm        
## 3         120 norm        
## 4         240 norm        
## 5          30 chi         
## 6          60 chi         
## 7         120 chi         
## 8         240 chi

Generate <- function(condition, fixed_objects = NULL) {
    N <- condition$sample_size
    dist <- condition$distribution
    if(dist == 'norm'){
        dat <- rnorm(N, mean = 3)
    } else if(dist == 'chi'){
        dat <- rchisq(N, df = 3)
    }
    dat
}

Define Bayes estimator of the mean with STAN

# STAN model
bmean_stan <- "
    data {
        int<lower=0> N;
        real x[N];
    }
    parameters {
        real mu;
        real<lower=0> sigma;
    }
    model {
        target += normal_lpdf(mu | 0, 10);  // weakly informative prior
        target += normal_lpdf(x | mu, sigma);
    }
"
# Save file
stan_path <- write_stan_file(bmean_stan)
mod <- cmdstan_model(stan_path)

## Warning in '/tmp/RtmpVU5PBw/model-7f13b1ada25b1.stan', line 4, column 8: Declaration
##     of arrays by placing brackets after a variable name is deprecated and
##     will be removed in Stan 2.32.0. Instead use the array keyword before the
##     type. This can be changed automatically using the auto-format flag to
##     stanc

Analyse <- function(condition, dat, fixed_objects = NULL) {
    mod <- fixed_objects$mod
    M0 <- mean(dat)
    M1 <- mean(dat, trim = .1)
    M2 <- mean(dat, trim = .2)
    med <- median(dat)
    stan_fit <- quiet(mod$sample(list(x = dat, N = length(dat)),
                                 refresh = 0, chains = 1, 
                                 show_messages = FALSE))
    MB <- stan_fit$summary("mu", mean)$mean[1]
    ret <- c(mean_no_trim = M0, mean_trim.1 = M1, 
             mean_trim.2 = M2, median = med, 
             bayes_mean = MB)
    ret
}

Summarise <- function(condition, results, fixed_objects = NULL) {
    obs_bias <- bias(results, parameter = 3)
    obs_RMSE <- RMSE(results, parameter = 3)
    ret <- c(bias = obs_bias, RMSE = obs_RMSE, RE = RE(obs_RMSE))
    ret
}

res <- runSimulation(Design, replications = 50, generate = Generate, 
                     analyse = Analyse, summarise = Summarise, 
                     parallel = TRUE,
                     ncores = min(2, parallel::detectCores()), 
                     fixed_objects = list(mod = mod),
                     packages = "cmdstanr")

## 
## Number of parallel clusters in use: 2

## 
## 
Design row: 1/8;   Started: Sun Jul 31 00:06:14 2022;   Total elapsed time: 0.00s 
## 
## 
Design row: 2/8;   Started: Sun Jul 31 00:06:21 2022;   Total elapsed time: 6.80s 
## 
## 
Design row: 3/8;   Started: Sun Jul 31 00:06:27 2022;   Total elapsed time: 12.65s 
## 
## 
Design row: 4/8;   Started: Sun Jul 31 00:06:34 2022;   Total elapsed time: 19.65s 
## 
## 
Design row: 5/8;   Started: Sun Jul 31 00:06:39 2022;   Total elapsed time: 25.48s 
## 
## 
Design row: 6/8;   Started: Sun Jul 31 00:06:45 2022;   Total elapsed time: 31.54s 
## 
## 
Design row: 7/8;   Started: Sun Jul 31 00:06:52 2022;   Total elapsed time: 38.26s 
## 
## 
Design row: 8/8;   Started: Sun Jul 31 00:06:59 2022;   Total elapsed time: 44.83s

## 
## Simulation complete. Total execution time: 50.90s

knitr::kable(res)