Merge pull request #9 from mtwesley/add-test-scripts

Add test scripts for custom censored distributions

tests/testthat/test-beta_censored.R

@@ -0,0 +1,45 @@
+# Test script for beta_censored distribution
+
+library(greta)
+library(testthat)
+
+test_that("beta_censored distribution works correctly", {
+  # Simulate data
+  set.seed(505)
+  n <- 100
+  true_alpha <- 2
+  true_beta <- 5
+  y <- rbeta(n, shape1 = true_alpha, shape2 = true_beta)
+
+  # Introduce interval censoring between 0.2 and 0.8
+  lower_bound <- 0.2
+  upper_bound <- 0.8
+  is_censored <- y > lower_bound & y < upper_bound
+  y_obs <- y
+  y_obs[is_censored] <- NA  # Interval censored data
+
+  # Data preparation
+  y_greta <- as_data(ifelse(is.na(y_obs), 0, y_obs))  # Placeholder for censored data
+  is_censored_greta <- as_data(as.numeric(is_censored))
+
+  # Define the model
+  alpha <- variable(lower = 0)
+  beta <- variable(lower = 0)
+
+  distribution(y_greta) <- beta_censored(
+    alpha = alpha,
+    beta = beta,
+    is_censored = is_censored_greta,
+    censoring_type = "interval",
+    lower = lower_bound,
+    upper = upper_bound,
+    dim = n
+  )
+
+  # Model fitting
+  m <- model(alpha, beta)
+  draws <- mcmc(m, n_samples = 1000)
+
+  # Output results
+  summary(draws)
+})

tests/testthat/test-exponential_censored.R

@@ -0,0 +1,40 @@
+# Test script for exponential_censored distribution
+
+library(greta)
+library(testthat)
+
+test_that("exponential_censored distribution works correctly", {
+  # Simulate data
+  set.seed(202)
+  n <- 100
+  true_rate <- 0.5
+  y <- rexp(n, rate = true_rate)
+
+  # Introduce left censoring at y < 0.5
+  censoring_threshold <- 0.5
+  is_censored <- y < censoring_threshold
+  y_obs <- ifelse(is_censored, censoring_threshold, y)
+
+  # Data preparation
+  y_greta <- as_data(y_obs)
+  is_censored_greta <- as_data(as.numeric(is_censored))
+
+  # Define the model
+  rate <- variable(lower = 0)
+
+  distribution(y_greta) <- exponential_censored(
+    rate = rate,
+    is_censored = is_censored_greta,
+    censoring_type = "left",
+    lower = NULL,
+    upper = NULL,
+    dim = n
+  )
+
+  # Model fitting
+  m <- model(rate)
+  draws <- mcmc(m, n_samples = 1000)
+
+  # Output results
+  summary(draws)
+})

tests/testthat/test-gamma_censored.R

@@ -0,0 +1,43 @@
+# Test script for gamma_censored distribution
+
+library(greta)
+library(testthat)
+
+test_that("gamma_censored distribution works correctly", {
+  # Simulate data
+  set.seed(101)
+  n <- 100
+  true_shape <- 2
+  true_rate <- 1
+  y <- rgamma(n, shape = true_shape, rate = true_rate)
+
+  # Introduce right censoring at y > 3
+  censoring_threshold <- 3
+  is_censored <- y > censoring_threshold
+  y_obs <- ifelse(is_censored, censoring_threshold, y)
+
+  # Data preparation
+  y_greta <- as_data(y_obs)
+  is_censored_greta <- as_data(as.numeric(is_censored))
+
+  # Define the model
+  shape <- variable(lower = 0)
+  rate <- variable(lower = 0)
+
+  distribution(y_greta) <- gamma_censored(
+    shape = shape,
+    rate = rate,
+    is_censored = is_censored_greta,
+    censoring_type = "right",
+    lower = NULL,
+    upper = NULL,
+    dim = n
+  )
+
+  # Model fitting
+  m <- model(shape, rate)
+  draws <- mcmc(m, n_samples = 1000)
+
+  # Output results
+  summary(draws)
+})

tests/testthat/test-lognormal_censored.R

@@ -0,0 +1,43 @@
+# Test script for lognormal_censored distribution
+
+library(greta)
+library(testthat)
+
+test_that("lognormal_censored distribution works correctly", {
+  # Simulate data
+  set.seed(456)
+  n <- 100
+  true_meanlog <- 0.5
+  true_sdlog <- 0.75
+  y <- rlnorm(n, meanlog = true_meanlog, sdlog = true_sdlog)
+
+  # Introduce left censoring at y < 1
+  censoring_threshold <- 1
+  is_censored <- y < censoring_threshold
+  y_obs <- ifelse(is_censored, censoring_threshold, y)
+
+  # Data preparation
+  y_greta <- as_data(y_obs)
+  is_censored_greta <- as_data(as.numeric(is_censored))
+
+  # Define the model
+  meanlog <- variable()
+  sdlog <- variable(lower = 0)
+
+  distribution(y_greta) <- lognormal_censored(
+    meanlog = meanlog,
+    sdlog = sdlog,
+    is_censored = is_censored_greta,
+    censoring_type = "left",
+    lower = NULL,
+    upper = NULL,
+    dim = n
+  )
+
+  # Model fitting
+  m <- model(meanlog, sdlog)
+  draws <- mcmc(m, n_samples = 1000)
+
+  # Output results
+  summary(draws)
+})

tests/testthat/test-normal_censored.R

@@ -0,0 +1,43 @@
+# Test script for normal_censored distribution
+
+library(greta)
+library(testthat)
+
+test_that("normal_censored distribution works correctly", {
+  # Simulate data
+  set.seed(123)
+  n <- 100
+  true_mean <- 2
+  true_sd <- 1
+  y <- rnorm(n, mean = true_mean, sd = true_sd)
+
+  # Introduce right censoring at y > 3
+  censoring_threshold <- 3
+  is_censored <- y > censoring_threshold
+  y_obs <- ifelse(is_censored, censoring_threshold, y)
+
+  # Data preparation
+  y_greta <- as_data(y_obs)
+  is_censored_greta <- as_data(as.numeric(is_censored))
+
+  # Define the model
+  mean <- variable()
+  sd <- variable(lower = 0)
+
+  distribution(y_greta) <- normal_censored(
+    mean = mean,
+    sd = sd,
+    is_censored = is_censored_greta,
+    censoring_type = "right",
+    lower = NULL,
+    upper = NULL,
+    dim = n
+  )
+
+  # Model fitting
+  m <- model(mean, sd)
+  draws <- mcmc(m, n_samples = 1000)
+
+  # Output results
+  summary(draws)
+})

tests/testthat/test-pareto_censored.R

@@ -0,0 +1,44 @@
+# Test script for pareto_censored distribution
+
+library(greta)
+library(testthat)
+
+test_that("pareto_censored distribution works correctly", {
+  # Simulate data
+  set.seed(404)
+  n <- 100
+  true_scale <- 1
+  true_alpha <- 2.5
+  library(VGAM)  # For rpareto
+  y <- rpareto(n, scale = true_scale, shape = true_alpha)
+
+  # Introduce left censoring at y < 2
+  censoring_threshold <- 2
+  is_censored <- y < censoring_threshold
+  y_obs <- ifelse(is_censored, censoring_threshold, y)
+
+  # Data preparation
+  y_greta <- as_data(y_obs)
+  is_censored_greta <- as_data(as.numeric(is_censored))
+
+  # Define the model
+  scale <- variable(lower = 0)
+  alpha <- variable(lower = 0)
+
+  distribution(y_greta) <- pareto_censored(
+    scale = scale,
+    alpha = alpha,
+    is_censored = is_censored_greta,
+    censoring_type = "left",
+    lower = NULL,
+    upper = NULL,
+    dim = n
+  )
+
+  # Model fitting
+  m <- model(scale, alpha)
+  draws <- mcmc(m, n_samples = 1000)
+
+  # Output results
+  summary(draws)
+})