## ----setup, include=FALSE-----------------------------------------------------
knitr::opts_chunk$set(echo = TRUE)
library(TKApprox)

## -----------------------------------------------------------------------------
# Define the exponential PDF
pdf_exp <- function(x, param) {
  dexp(x, rate = param)
}

# Define the exponential CDF
cdf_exp <- function(x, param) {
  pexp(x, rate = param)
}

## -----------------------------------------------------------------------------
# Gamma prior for the rate parameter
prior_spec <- list(
  rate = list(family = "gamma", hyperparameters = list(shape = 2, rate = 1))
)

## -----------------------------------------------------------------------------
set.seed(123)
data <- rexp(20, rate = 1.5)

## -----------------------------------------------------------------------------
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"
)

## -----------------------------------------------------------------------------
# Print basic results
print(fit)

# Detailed summary
summary(fit)

# Extract Bayes estimates
coef(fit)

# Extract covariance matrix
vcov(fit)

# Model comparison statistics
print_model_comparison(fit)

## -----------------------------------------------------------------------------
# Plot all diagnostics
plot(fit)

# Or select specific plots
plot(fit, which = 1)  # Posterior approximation
plot(fit, which = 2)  # Likelihood surface
plot(fit, which = 3)  # Prior vs posterior

## -----------------------------------------------------------------------------
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_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)
)
coef(fit_linex)

## -----------------------------------------------------------------------------
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)
)
coef(fit_gel)

## -----------------------------------------------------------------------------
# Create right-censored data
status <- c(1, 1, 0, 1, 0, 1, 1, 0, 1, 1)  # 1 = observed, 0 = censored

fit_censored <- 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
)

summary(fit_censored)

## -----------------------------------------------------------------------------
sensitivity <- tk_sensitivity(
  fit = fit,
  parameter_name = "rate",
  hyperparameter_name = "shape",
  hyperparameter_values = c(0.5, 1, 2, 5, 10)
)

print(sensitivity)
plot(sensitivity)

## -----------------------------------------------------------------------------
# Define 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])
}

# Independent priors
prior_spec_weibull <- list(
  shape = list(family = "gamma", hyperparameters = list(shape = 2, rate = 1)),
  scale = list(family = "gamma", hyperparameters = list(shape = 2, rate = 1))
)

# Generate Weibull data
set.seed(123)
data_weibull <- rweibull(20, shape = 2, scale = 1)

# Fit the model
fit_weibull <- tk_fit(
  data = data_weibull,
  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"
)

summary(fit_weibull)
plot(fit_weibull, which = 2)  # Likelihood surface for 2-parameter model

