Added stan_data to simulate_gastempt

R/simulate_gastempt.R

@@ -9,16 +9,20 @@
 #' @param student_t_df When NULL (default), Gaussian noise is added; when >= 2, Student_t distributed noise is added, which generates more realistic outliers. Values from 2 to 5 are useful, when higher values are used the result comes close to that of Gaussian noise. Values below 2 are truncated to 2.
 #' @param missing When 0 (default), all curves have the same number of data points. When > 0, this is the fraction of points that were removed randomly to simulate missing points. Maximum value is 0.5.
 #' @param model linexp(default) or powexp
+#' @param seed optional seed; not set if seed = NULL (default)
 #'
-#' @return A list with two elements:
+#' @return A list with 3 elements:
 #' \describe{
 #'   \item{record}{Data frame with columns
 #'     \code{record(chr), v0, tempt, kappa/beta} giving the effective
 #'     \code{linexp} or \code{powexp} parameters for the individual record.
 #'     \code{v0} is rounded to nearest integer.}
 #'   \item{data}{Data frame with columns
 #'     \code{record(chr), minute(dbl), vol(dbl)} giving the
-#'      time series and grouping parameters. \code{vol} is rounded to nearest integer.}
+#'      time series and grouping parameters. \code{vol} is rounded
+#'      to nearest integer.}
+#'   \item{stan_data}{A list for use as \code{data} in Stan-based fits
+#'   with elements \code{prior_v0, n, n_record, record, minute, volume}.}
 #'  }
 #'  A comment is attached to the return value that can be used as a title
 #' @examples
@@ -49,10 +53,13 @@ simulate_gastempt = function(
   noise = 20,
   student_t_df = NULL,
   missing = 0,
-  model = linexp) {
+  model = linexp,
+  seed = NULL) {
 
   # Hack to avoid notes
   vol = . = NULL
+  if (!is.null(seed))
+    set.seed(seed)
   # Only linexp and powexp are supported
   assert_that(identical(model, linexp) || identical(model, powexp))
   model_name = ifelse(identical(model, linexp), "linexp", "powexp")
@@ -132,7 +139,22 @@ simulate_gastempt = function(
           n_records, model_name, v0_mean, v0_std, tempt_mean, tempt_std, kb,
           sf, noise, 100*missing)
   comment(data) = comment
-  list(record = rec, data = data)
+  # Add list format for use with stan
+  #
+  data1 = data %>%
+    mutate(
+      record_i = as.integer(as.factor(data$record))
+    )
+
+  stan_data = list(
+    prior_v0 = median(data1$vol[data1$minute < 10]),
+    n = nrow(data1),
+    n_record = max(data1$record_i),
+    record = data1$record_i,
+    minute = data1$minute,
+    volume = data1$vol)
+
+  list(record = rec, data = data, stan_data = stan_data)
 }
 
 

tests/testthat/test-simulate_gastempt.R

@@ -6,7 +6,7 @@ library(assertthat)
 test_that("Default call of simulate_gastempt must return plausible values",{
   set.seed(4711)
   d = simulate_gastempt()
-  expect_equal(names(d), c("record", "data"))
+  expect_equal(names(d), c("record", "data", "stan_data"))
   expect_match(comment(d$data), "linexp")
   data = d$data
   record = d$record
@@ -20,7 +20,7 @@ test_that("Default call of simulate_gastempt must return plausible values",{
 test_that("Noise = 0 must issue a warning",{
   set.seed(4711)
   expect_warning(d <- simulate_gastempt(noise = 0), "zero")
-  expect_equal(names(d), c("record", "data"))
+  expect_equal(names(d), c("record", "data", "stan_data"))
   data = d$data
   record = d$record
   expect_equal(nrow(record), 10)
@@ -67,7 +67,7 @@ test_that(
   "Default call of simulate_gastempt for powexp must return plausible values",{
   set.seed(4711)
   d = simulate_gastempt(model = powexp)
-  expect_equal(names(d), c("record", "data"))
+  expect_equal(names(d), c("record", "data", "stan_data"))
   data = d$data
   record = d$record
   expect_equal(nrow(record), 10)
@@ -83,7 +83,7 @@ test_that("Noise = 0 must issue a warning and powexp should not overshoot",{
     noise = 0, model = powexp), "zero")
   data = d$data
   record = d$record
-  expect_equal(names(d), c("record", "data"))
+  expect_equal(names(d), c("record", "data", "stan_data"))
   expect_true(all(data$vol >= 0))
   expect_true(all(data$vol <= 400))
   # Without noise, first

tests/testthat/test-stan.R

@@ -19,17 +19,7 @@ test_that("stanmodels$linexp_gastro_1b exists", {
 })
 
 gastempt_data = function(){
-  set.seed(471)
-  d = simulate_gastempt(n_records = 6)$data
-  d$record_i =  as.integer(as.factor(d$record))
-
-  list(
-    prior_v0 = median(d$vol[d$minute < 10]),
-    n = nrow(d),
-    n_record = max(d$record_i),
-    record = d$record_i,
-    minute = d$minute,
-    volume = d$vol)
+  simulate_gastempt(n_records = 6, seed = 471)$stan_data
 }
 
 test_that("Direct use of sample model returns valid results", {
@@ -38,7 +28,7 @@ test_that("Direct use of sample model returns valid results", {
   stan_model = "../../exec/linexp_gastro_1b.stan"
   expect_true(file.exists(stan_model))
   rstan_options(auto_write = TRUE)
-  iter = 5000
+  iter = 500
   mr_b = stan(stan_model, data = data, chains = 4, iter = iter, seed = 4711,
               refresh = FALSE)
   expect_is(mr_b, "stanfit")