Simulation Studies

The hardware and bandwidth for this mirror is donated by METANET, the Webhosting and Full Service-Cloud Provider.
If you wish to report a bug, or if you are interested in having us mirror your free-software or open-source project, please feel free to contact us at mirror[@]metanet.ch.

Simulation Studies

Your Name

2026-07-31

Introduction

Simulation studies are essential for evaluating the performance of Bayesian estimation methods. This vignette demonstrates how to conduct simulation studies with TKApprox to assess bias, variance, mean squared error, and coverage probabilities of different estimators.

Basic Simulation Framework

A typical simulation study involves: 1. Generating data from a known distribution with known parameters 2. Fitting the model using TKApprox 3. Comparing estimates to true values 4. Repeating many times to assess performance

Example 1: Exponential Distribution with Complete Data

Setup

# Define exponential distribution
pdf_exp <- function(x, param) dexp(x, rate = param)
cdf_exp <- function(x, param) pexp(x, rate = param)

# Prior specification
prior_spec <- list(rate = list(family = "gamma", hyperparameters = list(shape = 2, rate = 1)))

# Simulation parameters
true_rate <- 1.5
sample_sizes <- c(20, 50, 100)
n_sim <- 10  # Number of simulations (small for fast vignette rendering)

Simulation Function

run_simulation <- function(n, true_rate, n_sim) {
  estimates_sel <- numeric(n_sim)
  estimates_linex <- numeric(n_sim)
  estimates_gel <- numeric(n_sim)
  
  for (i in 1:n_sim) {
    set.seed(i)
    data <- rexp(n, rate = true_rate)
    
    # SEL
    fit_sel <- tk_fit(
      data = data,
      censoring_scheme = "complete",
      pdf = pdf_exp,
      cdf = cdf_exp,
      prior_spec = prior_spec,
      initial_values = c(rate = 1),
      loss_function = "sel"
    )
    estimates_sel[i] <- coef(fit_sel)
    
    # LINEX
    fit_linex <- tk_fit(
      data = data,
      censoring_scheme = "complete",
      pdf = pdf_exp,
      cdf = cdf_exp,
      prior_spec = prior_spec,
      initial_values = c(rate = 1),
      loss_function = "linex",
      loss_params = list(c = 0.5)
    )
    estimates_linex[i] <- coef(fit_linex)
    
    # GEL
    fit_gel <- tk_fit(
      data = data,
      censoring_scheme = "complete",
      pdf = pdf_exp,
      cdf = cdf_exp,
      prior_spec = prior_spec,
      initial_values = c(rate = 1),
      loss_function = "gel",
      loss_params = list(q = 0.5)
    )
    estimates_gel[i] <- coef(fit_gel)
  }
  
  list(sel = estimates_sel, linex = estimates_linex, gel = estimates_gel)
}

Run Simulations

results <- lapply(sample_sizes, function(n) {
  run_simulation(n, true_rate, n_sim)
})
names(results) <- paste0("n_", sample_sizes)

Compute Performance Metrics

compute_metrics <- function(estimates, true_value) {
  bias <- mean(estimates) - true_value
  variance <- var(estimates)
  mse <- mean((estimates - true_value)^2)
  relative_bias <- bias / true_value
  rmse <- sqrt(mse)
  
  data.frame(
    bias = bias,
    variance = variance,
    mse = mse,
    relative_bias = relative_bias,
    rmse = rmse
  )
}

metrics <- lapply(results, function(res) {
  data.frame(
    Loss = c("SEL", "LINEX", "GEL"),
    rbind(
      compute_metrics(res$sel, true_rate),
      compute_metrics(res$linex, true_rate),
      compute_metrics(res$gel, true_rate)
    )
  )
})

Display Results

for (i in seq_along(sample_sizes)) {
  cat("\n=== Sample Size:", sample_sizes[i], "===\n")
  print(metrics[[i]])
}

## 
## === Sample Size: 20 ===
##    Loss      bias   variance       mse relative_bias      rmse
## 1   SEL 0.1626120 0.10584171 0.1217002    0.10840797 0.3488555
## 2 LINEX 0.1306818 0.09732687 0.1046719    0.08712118 0.3235304
## 3   GEL 0.1060642 0.09876448 0.1001376    0.07070944 0.3164453
## 
## === Sample Size: 50 ===
##    Loss       bias   variance        mse relative_bias      rmse
## 1   SEL 0.06754210 0.03417966 0.03532363    0.04502807 0.1879458
## 2 LINEX 0.05565425 0.03309691 0.03288462    0.03710283 0.1813412
## 3   GEL 0.04495315 0.03320166 0.03190228    0.02996876 0.1786121
## 
## === Sample Size: 100 ===
##    Loss       bias    variance         mse relative_bias       rmse
## 1   SEL 0.03033529 0.005712953 0.006061887    0.02022353 0.07785812
## 2 LINEX 0.02459912 0.005627390 0.005669767    0.01639941 0.07529786
## 3   GEL 0.01908762 0.005629284 0.005430693    0.01272508 0.07369324

Visualize Results

# Plot bias vs sample size
bias_sel <- sapply(metrics, function(m) m$bias[1])
bias_linex <- sapply(metrics, function(m) m$bias[2])
bias_gel <- sapply(metrics, function(m) m$bias[3])

plot(sample_sizes, bias_sel, type = "b", pch = 19, col = "blue",
     ylim = range(c(bias_sel, bias_linex, bias_gel)),
     xlab = "Sample Size", ylab = "Bias",
     main = "Bias vs Sample Size")
lines(sample_sizes, bias_linex, type = "b", pch = 19, col = "red")
lines(sample_sizes, bias_gel, type = "b", pch = 19, col = "green")
legend("topright", legend = c("SEL", "LINEX", "GEL"),
       col = c("blue", "red", "green"), pch = 19, lty = 1)

# Plot MSE vs sample size
mse_sel <- sapply(metrics, function(m) m$mse[1])
mse_linex <- sapply(metrics, function(m) m$mse[2])
mse_gel <- sapply(metrics, function(m) m$mse[3])

plot(sample_sizes, mse_sel, type = "b", pch = 19, col = "blue",
     ylim = range(c(mse_sel, mse_linex, mse_gel)),
     xlab = "Sample Size", ylab = "MSE",
     main = "MSE vs Sample Size")
lines(sample_sizes, mse_linex, type = "b", pch = 19, col = "red")
lines(sample_sizes, mse_gel, type = "b", pch = 19, col = "green")
legend("topright", legend = c("SEL", "LINEX", "GEL"),
       col = c("blue", "red", "green"), pch = 19, lty = 1)

Example 2: Censoring Schemes Comparison

Setup

# Simulation parameters
true_rate <- 1.5
n <- 50
n_sim <- 10

# Censoring proportions
censoring_props <- c(0.2, 0.4, 0.6)

Simulation Function for Right Censoring

run_censoring_simulation <- function(censoring_prop, n, true_rate, n_sim) {
  estimates <- numeric(n_sim)
  
  for (i in 1:n_sim) {
    set.seed(i)
    data <- rexp(n, rate = true_rate)
    
    # Apply right censoring
    censoring_time <- quantile(data, 1 - censoring_prop)
    status <- as.numeric(data <= censoring_time)
    
    fit <- tk_fit(
      data = data,
      censoring_scheme = "right-censored",
      pdf = pdf_exp,
      cdf = cdf_exp,
      prior_spec = prior_spec,
      initial_values = c(rate = 1),
      loss_function = "sel",
      status = status
    )
    
    estimates[i] <- coef(fit)
  }
  
  estimates
}

Run Simulations

censoring_results <- lapply(censoring_props, function(prop) {
  run_censoring_simulation(prop, n, true_rate, n_sim)
})
names(censoring_results) <- paste0("censoring_", censoring_props)

Compute Metrics

censoring_metrics <- lapply(censoring_results, function(est) {
  compute_metrics(est, true_rate)
})

Display Results

censoring_df <- do.call(rbind, censoring_metrics)
censoring_df$censoring_prop <- censoring_props
print(censoring_df)

##                     bias    variance        mse relative_bias      rmse
## censoring_0.2 -0.2338869 0.022298444 0.07477166    -0.1559246 0.2734441
## censoring_0.4 -0.5353080 0.012945138 0.29820529    -0.3568720 0.5460818
## censoring_0.6 -0.8367103 0.006119814 0.70559202    -0.5578069 0.8399952
##               censoring_prop
## censoring_0.2            0.2
## censoring_0.4            0.4
## censoring_0.6            0.6

Visualize

plot(censoring_props, censoring_metrics$bias, type = "b", pch = 19,
     xlab = "Censoring Proportion", ylab = "Bias",
     main = "Bias vs Censoring Proportion")
abline(h = 0, col = "red", lty = 2)

plot(censoring_props, censoring_metrics$mse, type = "b", pch = 19,
     xlab = "Censoring Proportion", ylab = "MSE",
     main = "MSE vs Censoring Proportion")

Example 3: Prior Sensitivity Simulation

Setup

# Different prior specifications
prior_specs <- list(
  weak = list(rate = list(family = "gamma", hyperparameters = list(shape = 0.1, rate = 0.1))),
  moderate = list(rate = list(family = "gamma", hyperparameters = list(shape = 2, rate = 1))),
  strong = list(rate = list(family = "gamma", hyperparameters = list(shape = 10, rate = 5)))
)

true_rate <- 1.5
n <- 50
n_sim <- 10

Simulation Function

run_prior_simulation <- function(prior_spec, n, true_rate, n_sim) {
  estimates <- numeric(n_sim)
  
  for (i in 1:n_sim) {
    set.seed(i)
    data <- rexp(n, rate = true_rate)
    
    fit <- tk_fit(
      data = data,
      censoring_scheme = "complete",
      pdf = pdf_exp,
      cdf = cdf_exp,
      prior_spec = prior_spec,
      initial_values = c(rate = 1),
      loss_function = "sel"
    )
    
    estimates[i] <- coef(fit)
  }
  
  estimates
}

Run Simulations

prior_results <- lapply(prior_specs, function(pspec) {
  run_prior_simulation(pspec, n, true_rate, n_sim)
})

Compute Metrics

prior_metrics <- lapply(prior_results, function(est) {
  compute_metrics(est, true_rate)
})

Display Results

prior_df <- do.call(rbind, prior_metrics)
prior_df$prior_strength <- names(prior_specs)
print(prior_df)

##                bias   variance        mse relative_bias      rmse
## weak     0.05294538 0.03571617 0.03494776    0.03529692 0.1869432
## moderate 0.06754210 0.03417966 0.03532363    0.04502807 0.1879458
## strong   0.11215924 0.02797773 0.03775965    0.07477282 0.1943184
##          prior_strength
## weak               weak
## moderate       moderate
## strong           strong

Example 4: Coverage Probability

Setup

true_rate <- 1.5
n <- 50
n_sim <- 20
alpha <- 0.05  # For 95% credible intervals

Simulation Function

run_coverage_simulation <- function(n, true_rate, n_sim, alpha) {
  coverage_count <- 0
  
  for (i in 1:n_sim) {
    set.seed(i)
    data <- rexp(n, rate = true_rate)
    
    fit <- tk_fit(
      data = data,
      censoring_scheme = "complete",
      pdf = pdf_exp,
      cdf = cdf_exp,
      prior_spec = prior_spec,
      initial_values = c(rate = 1),
      loss_function = "sel"
    )
    
    ci <- fit$credible_intervals
    lower <- ci[1, 1]
    upper <- ci[1, 2]
    
    if (true_rate >= lower && true_rate <= upper) {
      coverage_count <- coverage_count + 1
    }
  }
  
  coverage_count / n_sim
}

Run Simulation

coverage_prob <- run_coverage_simulation(n, true_rate, n_sim, alpha)
cat("Coverage probability:", coverage_prob, "\n")

## Coverage probability: 1

cat("Nominal coverage:", 1 - alpha, "\n")

## Nominal coverage: 0.95

Example 5: Multi-Parameter Model Simulation

Setup

# Weibull distribution
pdf_weibull <- function(x, param) dweibull(x, shape = param[1], scale = param[2])
cdf_weibull <- function(x, param) pweibull(x, shape = param[1], scale = param[2])

prior_spec_weibull <- list(
  shape = list(family = "gamma", hyperparameters = list(shape = 2, rate = 1)),
  scale = list(family = "gamma", hyperparameters = list(shape = 2, rate = 1))
)

true_shape <- 2
true_scale <- 1
n <- 50
n_sim <- 10

Simulation Function

run_weibull_simulation <- function(n, true_shape, true_scale, n_sim) {
  shape_estimates <- numeric(n_sim)
  scale_estimates <- numeric(n_sim)
  
  for (i in 1:n_sim) {
    set.seed(i)
    data <- rweibull(n, shape = true_shape, scale = true_scale)
    
    fit <- tk_fit(
      data = data,
      censoring_scheme = "complete",
      pdf = pdf_weibull,
      cdf = cdf_weibull,
      prior_spec = prior_spec_weibull,
      initial_values = c(shape = 1.5, scale = 1),
      loss_function = "sel"
    )
    
    shape_estimates[i] <- coef(fit)[1]
    scale_estimates[i] <- coef(fit)[2]
  }
  
  list(shape = shape_estimates, scale = scale_estimates)
}

Run Simulation

weibull_results <- run_weibull_simulation(n, true_shape, true_scale, n_sim)

Compute Metrics

shape_metrics <- compute_metrics(weibull_results$shape, true_shape)
scale_metrics <- compute_metrics(weibull_results$scale, true_scale)

cat("Shape parameter:\n")

## Shape parameter:

print(shape_metrics)

##          bias   variance        mse relative_bias      rmse
## 1 -0.02605768 0.09525672 0.08641005   -0.01302884 0.2939559

cat("\nScale parameter:\n")

## 
## Scale parameter:

print(scale_metrics)

##          bias    variance         mse relative_bias       rmse
## 1 -0.01898825 0.004473991 0.004387145   -0.01898825 0.06623553

Parallel Simulation

For large simulation studies, you can use parallel processing:

# Note: This requires the parallel package
library(parallel)

run_parallel_simulation <- function(n, true_rate, n_sim, n_cores = 4) {
  cl <- makeCluster(n_cores)
  
  results <- parLapply(cl, 1:n_sim, function(i) {
    set.seed(i)
    data <- rexp(n, rate = true_rate)
    
    fit <- tk_fit(
      data = data,
      censoring_scheme = "complete",
      pdf = pdf_exp,
      cdf = cdf_exp,
      prior_spec = prior_spec,
      initial_values = c(rate = 1),
      loss_function = "sel"
    )
    
    coef(fit)
  })
  
  stopCluster(cl)
  unlist(results)
}

Tips for Simulation Studies

Set seeds: Always set seeds for reproducibility
Number of simulations: Use at least 100-1000 simulations for stable estimates
Sample sizes: Test a range of sample sizes to assess asymptotic behavior
Convergence: Monitor optimization convergence in simulations
Parallel processing: Use parallel processing for large simulation studies
Store results: Save simulation results for later analysis
Visualization: Always visualize simulation results

Interpreting Simulation Results

Bias

Small bias: Estimates are approximately unbiased
Large bias: May indicate prior influence or model misspecification
Direction of bias: Systematic over- or under-estimation

Variance

Low variance: Precise estimates
High variance: May need larger sample sizes or stronger priors

MSE

Combines bias and variance: MSE = Variance + Bias²
Lower MSE indicates better overall performance

Coverage Probability

Close to nominal: Credible intervals are well-calibrated
Below nominal: Intervals are too narrow (anti-conservative)
Above nominal: Intervals are too wide (conservative)

Next Steps

See “Real Data Example” for application to actual data
See “Loss Functions” for comparing estimators in simulation

These binaries (installable software) and packages are in development.
They may not be fully stable and should be used with caution. We make no claims about them.