---
title: "Simulation Studies"
author: "Your Name"
date: "`r Sys.Date()`"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Simulation Studies}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
library(TKApprox)
```

## 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

```{r}
# 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

```{r}
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

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

### Compute Performance Metrics

```{r}
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

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

### Visualize Results

```{r}
# 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

```{r}
# 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

```{r}
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

```{r}
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

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

### Display Results

```{r}
censoring_df <- do.call(rbind, censoring_metrics)
censoring_df$censoring_prop <- censoring_props
print(censoring_df)
```

### Visualize

```{r}
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

```{r}
# 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

```{r}
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

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

### Compute Metrics

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

### Display Results

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

## Example 4: Coverage Probability

### Setup

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

### Simulation Function

```{r}
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

```{r}
coverage_prob <- run_coverage_simulation(n, true_rate, n_sim, alpha)
cat("Coverage probability:", coverage_prob, "\n")
cat("Nominal coverage:", 1 - alpha, "\n")
```

## Example 5: Multi-Parameter Model Simulation

### Setup

```{r}
# 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

```{r}
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

```{r}
weibull_results <- run_weibull_simulation(n, true_shape, true_scale, n_sim)
```

### Compute Metrics

```{r}
shape_metrics <- compute_metrics(weibull_results$shape, true_shape)
scale_metrics <- compute_metrics(weibull_results$scale, true_scale)

cat("Shape parameter:\n")
print(shape_metrics)
cat("\nScale parameter:\n")
print(scale_metrics)
```

## Parallel Simulation

For large simulation studies, you can use parallel processing:

```{r}
# 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

1. **Set seeds**: Always set seeds for reproducibility
2. **Number of simulations**: Use at least 100-1000 simulations for stable estimates
3. **Sample sizes**: Test a range of sample sizes to assess asymptotic behavior
4. **Convergence**: Monitor optimization convergence in simulations
5. **Parallel processing**: Use parallel processing for large simulation studies
6. **Store results**: Save simulation results for later analysis
7. **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