## ----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_spec <- list(rate = list(family = "gamma", hyperparameters = list(shape = 2, rate = 1)))

set.seed(123)
data <- rexp(20, rate = 1.5)

# Fit with SEL (default)
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"
)

coef(fit_sel)

## -----------------------------------------------------------------------------
# Fit with LINEX loss (c = 0.5, penalizes overestimation)
fit_linex_pos <- 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)
)

# Fit with LINEX loss (c = -0.5, penalizes underestimation)
fit_linex_neg <- 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)
)

# Compare estimates
data.frame(
  SEL = coef(fit_sel),
  LINEX_c_0.5 = coef(fit_linex_pos),
  LINEX_c_neg0.5 = coef(fit_linex_neg)
)

## -----------------------------------------------------------------------------
# Examine sensitivity to LINEX parameter c
c_values <- c(-2, -1, -0.5, -0.1, 0.1, 0.5, 1, 2)
linex_estimates <- sapply(c_values, function(c) {
  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 = "linex",
    loss_params = list(c = c)
  )
  coef(fit)
})

plot(c_values, linex_estimates, type = "b", pch = 19,
     xlab = "LINEX parameter c", ylab = "Estimate",
     main = "LINEX Estimates vs c")
abline(h = coef(fit_sel), col = "red", lty = 2)
legend("topright", legend = c("LINEX", "SEL"), col = c("black", "red"),
       pch = c(19, NA), lty = c(1, 2))

## -----------------------------------------------------------------------------
# Fit with GEL (q = 0.5)
fit_gel_pos <- 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)
)

# Fit with GEL (q = -0.5)
fit_gel_neg <- 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)
)

# Compare estimates
data.frame(
  SEL = coef(fit_sel),
  GEL_q_0.5 = coef(fit_gel_pos),
  GEL_q_neg0.5 = coef(fit_gel_neg)
)

## -----------------------------------------------------------------------------
# Fit with precautionary loss
fit_precautionary <- tk_fit(
  data = data,
  censoring_scheme = "complete",
  pdf = pdf_exp,
  cdf = cdf_exp,
  prior_spec = prior_spec,
  initial_values = c(rate = 1),
  loss_function = "precautionary"
)

# Compare with SEL
data.frame(
  SEL = coef(fit_sel),
  Precautionary = coef(fit_precautionary)
)

## -----------------------------------------------------------------------------
# Fit with weighted SEL
fit_weighted <- tk_fit(
  data = data,
  censoring_scheme = "complete",
  pdf = pdf_exp,
  cdf = cdf_exp,
  prior_spec = prior_spec,
  initial_values = c(rate = 1),
  loss_function = "weighted-sel"
)

# Compare with SEL
data.frame(
  SEL = coef(fit_sel),
  Weighted_SEL = coef(fit_weighted)
)

## -----------------------------------------------------------------------------
# Example: Estimate the median of the posterior
# For exponential distribution, median = log(2)/rate
# We want to estimate log(rate) instead of rate directly
custom_g <- function(param) log(param[1])

fit_custom <- tk_fit(
  data = data,
  censoring_scheme = "complete",
  pdf = pdf_exp,
  cdf = cdf_exp,
  prior_spec = prior_spec,
  initial_values = c(rate = 1),
  loss_function = "custom",
  custom_g = custom_g
)

# The estimate is E[log(rate) | x]
# Transform back to rate scale
rate_estimate <- exp(coef(fit_custom))

data.frame(
  SEL_rate = coef(fit_sel),
  Custom_log_rate = rate_estimate
)

## -----------------------------------------------------------------------------
# Fit with all loss functions
fit_sel <- tk_fit(data, "complete", pdf_exp, cdf_exp, prior_spec,
                  initial_values = c(rate = 1), loss_function = "sel")

fit_linex <- tk_fit(data, "complete", pdf_exp, cdf_exp, prior_spec,
                    initial_values = c(rate = 1), loss_function = "linex", loss_params = list(c = 0.5))

fit_gel <- tk_fit(data, "complete", pdf_exp, cdf_exp, prior_spec,
                  initial_values = c(rate = 1), loss_function = "gel", loss_params = list(q = 0.5))

fit_precautionary <- tk_fit(data, "complete", pdf_exp, cdf_exp, prior_spec,
                            initial_values = c(rate = 1), loss_function = "precautionary")

fit_weighted <- tk_fit(data, "complete", pdf_exp, cdf_exp, prior_spec,
                       initial_values = c(rate = 1), loss_function = "weighted-sel")

# Compare all estimates
comparison <- data.frame(
  Loss_Function = c("SEL", "LINEX (c=0.5)", "GEL (q=0.5)", "Precautionary", "Weighted SEL"),
  Estimate = c(coef(fit_sel), coef(fit_linex), coef(fit_gel),
               coef(fit_precautionary), coef(fit_weighted))
)

print(comparison)

# Visual comparison
barplot(comparison$Estimate, names.arg = comparison$Loss_Function,
        main = "Bayes Estimates Under Different Loss Functions",
        ylab = "Rate Estimate", col = "steelblue")
abline(h = 1.5, col = "red", lty = 2)  # True value
legend("topright", legend = "True value", col = "red", lty = 2)

## -----------------------------------------------------------------------------
# Two-parameter Weibull example
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 <- list(
  shape = list(family = "gamma", hyperparameters = list(shape = 2, rate = 1)),
  scale = list(family = "gamma", hyperparameters = list(shape = 2, rate = 1))
)

set.seed(123)
data_weibull <- rweibull(20, shape = 2, scale = 1)

# Fit with SEL (applies to both parameters)
fit_weibull <- tk_fit(
  data = data_weibull,
  censoring_scheme = "complete",
  pdf = pdf_weibull,
  cdf = cdf_weibull,
  prior_spec = prior_spec,
  initial_values = c(shape = 1.5, scale = 1),
  loss_function = "sel"
)

coef(fit_weibull)