## ----setup, include=FALSE-----------------------------------------------------
knitr::opts_chunk$set(echo = TRUE)
library(TKApprox)

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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