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)
sample_size distribution bias.mean_no_trim bias.mean_trim.1 bias.mean_trim.2 bias.median bias.bayes_mean RMSE.mean_no_trim RMSE.mean_trim.1 RMSE.mean_trim.2 RMSE.median RMSE.bayes_mean RE.mean_no_trim RE.mean_trim.1 RE.mean_trim.2 RE.median RE.bayes_mean REPLICATIONS SIM_TIME COMPLETED SEED
30 norm -0.0118190 0.0014862 0.0126367 0.0417798 -0.0124052 0.1925357 0.1962211 0.1974239 0.2340095 0.1931643 1 1.038649 1.051421 1.477217 1.0065399 50 6.802 Sun Jul 31 00:06:21 2022 919770021
60 norm 0.0138500 0.0123798 0.0159013 0.0143448 0.0137575 0.1294069 0.1322249 0.1340543 0.1620916 0.1296239 1 1.044027 1.073117 1.568940 1.0033577 50 5.844 Sun Jul 31 00:06:27 2022 1147125301
120 norm 0.0248642 0.0251620 0.0262231 0.0309249 0.0243947 0.0903066 0.0925141 0.0966194 0.1155442 0.0896688 1 1.049486 1.144695 1.637034 0.9859249 50 7.003 Sun Jul 31 00:06:34 2022 165936352
240 norm -0.0124981 -0.0135860 -0.0141574 -0.0148045 -0.0123418 0.0681080 0.0685800 0.0715182 0.0827685 0.0688086 1 1.013908 1.102646 1.476840 1.0206777 50 5.828 Sun Jul 31 00:06:39 2022 1169202745
30 chi 0.0300504 -0.2925417 -0.4179901 -0.4775784 0.0219484 0.4460887 0.5256165 0.6078532 0.7030688 0.4428674 1 1.388339 1.856757 2.484010 0.9856100 50 6.064 Sun Jul 31 00:06:45 2022 666552130
60 chi 0.0302642 -0.3434107 -0.5061163 -0.6635279 0.0306421 0.3189606 0.4417358 0.5724486 0.7151231 0.3217270 1 1.918011 3.221060 5.026751 1.0174211 50 6.715 Sun Jul 31 00:06:52 2022 245891705
120 chi 0.0297133 -0.3006558 -0.4464628 -0.5809295 0.0288748 0.2420612 0.3748495 0.4983842 0.6294695 0.2401434 1 2.398079 4.239144 6.762371 0.9842171 50 6.571 Sun Jul 31 00:06:59 2022 1719562701
240 chi 0.0170796 -0.3419476 -0.4842108 -0.6231752 0.0166940 0.1400220 0.3683362 0.5032654 0.6455598 0.1398539 1 6.919842 12.918189 21.255941 0.9976008 50 6.076 Sun Jul 31 00:07:05 2022 1548188555
Hok Chio (Mark) Lai 黎學昭
Hok Chio (Mark) Lai 黎學昭
Assistant Professor of Psychology (Quantitative Methods)

My research interests include statistics, multilevel and latent variable models, and psychometrics.

comments powered by Disqus

Related