# Using cmdstanr in SimDesign

Programming
Statistics
Author

Mark Lai

Published

July 14, 2021

Modified

November 3, 2023

``````library(SimDesign)
library(cmdstanr)``````

[Update: Use parallel computing with two cores.]

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;
array[N] real x;
}
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)``````
``````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;   RAM used: 68.2 Mb;   Total elapsed time: 0.00s
Conditions: sample_size=30, distribution=norm

Design row: 2/8;   RAM used: 69 Mb;   Total elapsed time: 5.22s
Conditions: sample_size=60, distribution=norm

Design row: 3/8;   RAM used: 69 Mb;   Total elapsed time: 10.26s
Conditions: sample_size=120, distribution=norm

Design row: 4/8;   RAM used: 69 Mb;   Total elapsed time: 15.12s
Conditions: sample_size=240, distribution=norm

Design row: 5/8;   RAM used: 69 Mb;   Total elapsed time: 20.10s
Conditions: sample_size=30, distribution=chi

Design row: 6/8;   RAM used: 69 Mb;   Total elapsed time: 25.06s
Conditions: sample_size=60, distribution=chi

Design row: 7/8;   RAM used: 69 Mb;   Total elapsed time: 29.79s
Conditions: sample_size=120, distribution=chi

Design row: 8/8;   RAM used: 69 Mb;   Total elapsed time: 34.88s
Conditions: sample_size=240, distribution=chi
``````
``````
Simulation complete. Total execution time: 39.69s``````
``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 RAM_USED SEED COMPLETED
30 norm 0.0137592 0.0072189 0.0051610 -0.0088380 0.0134503 0.1831414 0.2016381 0.2068831 0.2458775 0.1837847 1 1.212194 1.276078 1.802455 1.0070379 50 5.223 69 Mb 77254733 Thu Mar 21 10:32:12 2024
60 norm -0.0161554 -0.0167838 -0.0165226 -0.0155586 -0.0181256 0.1171312 0.1192122 0.1279795 0.1617023 0.1176554 1 1.035847 1.193811 1.905843 1.0089697 50 5.041 69 Mb 745175193 Thu Mar 21 10:32:17 2024
120 norm -0.0135750 -0.0129037 -0.0129281 -0.0094663 -0.0135980 0.0763685 0.0843901 0.0888738 0.1060062 0.0766066 1 1.221109 1.354312 1.926789 1.0062448 50 4.852 69 Mb 796799714 Thu Mar 21 10:32:22 2024
240 norm 0.0112872 0.0172670 0.0194530 0.0195026 0.0113646 0.0661500 0.0699354 0.0717609 0.0786083 0.0664199 1 1.117724 1.176837 1.412139 1.0081765 50 4.982 69 Mb 1458160184 Thu Mar 21 10:32:27 2024
30 chi 0.0141920 -0.3026217 -0.4385604 -0.5725934 0.0047972 0.4931700 0.5686208 0.6421671 0.7338760 0.4913024 1 1.329389 1.695520 2.214380 0.9924404 50 4.962 69 Mb 153258056 Thu Mar 21 10:32:32 2024
60 chi 0.0219076 -0.3127238 -0.4623196 -0.6273301 0.0210254 0.2738739 0.4046596 0.5340195 0.6891814 0.2767046 1 2.183125 3.802007 6.332363 1.0207785 50 4.729 69 Mb 838557643 Thu Mar 21 10:32:37 2024
120 chi -0.0170950 -0.3392576 -0.4702292 -0.5832027 -0.0186802 0.2179343 0.3971243 0.5138743 0.6278862 0.2178307 1 3.320487 5.559843 8.300623 0.9990498 50 5.093 69 Mb 1973331824 Thu Mar 21 10:32:42 2024
240 chi 0.0242644 -0.3314124 -0.4775029 -0.6166151 0.0238784 0.1424573 0.3627346 0.5014147 0.6403622 0.1439531 1 6.483486 12.388673 20.206075 1.0211110 50 4.808 69 Mb 1726109243 Thu Mar 21 10:32:47 2024