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

## -----------------------------------------------------------------------------
# Gamma prior for exponential rate parameter
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 <- 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"
)

summary(fit)

## -----------------------------------------------------------------------------
# Normal prior for log-normal meanlog parameter
pdf_lognormal <- function(x, param) dlnorm(x, meanlog = param[1], sdlog = param[2])
cdf_lognormal <- function(x, param) plnorm(x, meanlog = param[1], sdlog = param[2])

prior_spec <- list(
  meanlog = list(family = "normal", hyperparameters = list(mean = 0, sd = 1)),
  sdlog = list(family = "gamma", hyperparameters = list(shape = 2, rate = 1))
)

set.seed(123)
data <- rlnorm(20, meanlog = 0, sdlog = 0.5)

fit <- tk_fit(
  data = data,
  censoring_scheme = "complete",
  pdf = pdf_lognormal,
  cdf = cdf_lognormal,
  prior_spec = prior_spec,
  initial_values = c(meanlog = 0, sdlog = 0.5),
  loss_function = "sel"
)

summary(fit)

## -----------------------------------------------------------------------------
# Beta prior for probability parameter
pdf_bernoulli <- function(x, param) {
  p <- param[1]
  ifelse(x == 1, p, 1 - p)
}

cdf_bernoulli <- function(x, param) {
  p <- param[1]
  ifelse(x == 0, 1 - p, 1)
}

prior_spec <- list(
  p = list(family = "beta", hyperparameters = list(shape1 = 2, shape2 = 2))
)

# Bernoulli data
set.seed(123)
data <- rbinom(20, size = 1, prob = 0.6)

fit <- tk_fit(
  data = data,
  censoring_scheme = "complete",
  pdf = pdf_bernoulli,
  cdf = cdf_bernoulli,
  prior_spec = prior_spec,
  initial_values = c(p = 0.5),
  loss_function = "sel"
)

summary(fit)

## -----------------------------------------------------------------------------
# Uniform prior for Weibull shape parameter
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 = "uniform", hyperparameters = list(lower = 0.1, upper = 10)),
  scale = list(family = "gamma", hyperparameters = list(shape = 2, rate = 1))
)

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

fit <- tk_fit(
  data = data,
  censoring_scheme = "complete",
  pdf = pdf_weibull,
  cdf = cdf_weibull,
  prior_spec = prior_spec,
  initial_values = c(shape = 1.5, scale = 1),
  loss_function = "sel"
)

summary(fit)

## -----------------------------------------------------------------------------
# Exponential prior for Poisson rate
pdf_poisson <- function(x, param) dpois(x, lambda = param[1])
cdf_poisson <- function(x, param) ppois(x, lambda = param[1])

prior_spec <- list(
  lambda = list(family = "exponential", hyperparameters = list(rate = 1))
)

set.seed(123)
data <- rpois(20, lambda = 3)

fit <- tk_fit(
  data = data,
  censoring_scheme = "complete",
  pdf = pdf_poisson,
  cdf = cdf_poisson,
  prior_spec = prior_spec,
  initial_values = c(lambda = 2),
  loss_function = "sel"
)

summary(fit)

## -----------------------------------------------------------------------------
# Log-Normal prior for Pareto scale parameter
pdf_pareto <- function(x, param) {
  xm <- param[1]
  alpha <- param[2]
  ifelse(x >= xm, (alpha * xm^alpha) / (x^(alpha + 1)), 0)
}

cdf_pareto <- function(x, param) {
  xm <- param[1]
  alpha <- param[2]
  ifelse(x >= xm, 1 - (xm / x)^alpha, 0)
}

prior_spec <- list(
  xm = list(family = "lognormal", hyperparameters = list(meanlog = 0, sdlog = 0.5)),
  alpha = list(family = "gamma", hyperparameters = list(shape = 2, rate = 1))
)

set.seed(123)
data <- (1 / (1 - runif(20)))^(1/2)  # Pareto(1, 2)

fit <- tk_fit(
  data = data,
  censoring_scheme = "complete",
  pdf = pdf_pareto,
  cdf = cdf_pareto,
  prior_spec = prior_spec,
  initial_values = c(xm = 0.5, alpha = 1.5),
  loss_function = "sel"
)

summary(fit)

## -----------------------------------------------------------------------------
# Weibull prior for gamma shape parameter
pdf_gamma <- function(x, param) dgamma(x, shape = param[1], rate = param[2])
cdf_gamma <- function(x, param) pgamma(x, shape = param[1], rate = param[2])

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

set.seed(123)
data <- rgamma(20, shape = 2, rate = 1.5)

fit <- tk_fit(
  data = data,
  censoring_scheme = "complete",
  pdf = pdf_gamma,
  cdf = cdf_gamma,
  prior_spec = prior_spec,
  initial_values = c(shape = 1.5, rate = 1),
  loss_function = "sel"
)

summary(fit)

## -----------------------------------------------------------------------------
# Inverse Gamma prior for normal variance
pdf_normal <- function(x, param) dnorm(x, mean = param[1], sd = sqrt(param[2]))
cdf_normal <- function(x, param) pnorm(x, mean = param[1], sd = sqrt(param[2]))

prior_spec <- list(
  mean = list(family = "normal", hyperparameters = list(mean = 0, sd = 10)),
  variance = list(family = "invgamma", hyperparameters = list(shape = 2, scale = 1))
)

set.seed(123)
data <- rnorm(20, mean = 0, sd = 2)

fit <- tk_fit(
  data = data,
  censoring_scheme = "complete",
  pdf = pdf_normal,
  cdf = cdf_normal,
  prior_spec = prior_spec,
  initial_values = c(mean = 0, variance = 4),
  loss_function = "sel"
)

summary(fit)

## -----------------------------------------------------------------------------
# Two-parameter 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 for each parameter
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 <- rweibull(20, shape = 2, scale = 1)

fit <- tk_fit(
  data = data,
  censoring_scheme = "complete",
  pdf = pdf_weibull,
  cdf = cdf_weibull,
  prior_spec = prior_spec,
  initial_values = c(shape = 1.5, scale = 1),
  loss_function = "sel"
)

summary(fit)

## -----------------------------------------------------------------------------
# Custom prior function
custom_logprior <- function(param) {
  # Example: hierarchical prior
  # param[1] = theta, param[2] = hyperparameter
  theta <- param[1]
  hyper <- param[2]
  
  # Prior for theta given hyper
  log_prior_theta <- dnorm(theta, mean = 0, sd = hyper, log = TRUE)
  
  # Prior for hyper
  log_prior_hyper <- dgamma(hyper, shape = 2, rate = 1, log = TRUE)
  
  log_prior_theta + log_prior_hyper
}

# Use custom prior in tk_fit
fit <- tk_fit(
  data = data,
  censoring_scheme = "complete",
  pdf = pdf_exp,
  cdf = cdf_exp,
  prior_spec = custom_logprior,
  initial_values = c(rate = 1),
  loss_function = "sel"
)

## -----------------------------------------------------------------------------
fit_flat <- tk_fit(
  data = data,
  censoring_scheme = "complete",
  pdf = pdf_exp,
  cdf = cdf_exp,
  prior_spec = NULL,  # Flat prior
  initial_values = c(rate = 1),
  loss_function = "sel"
)

summary(fit_flat)

## -----------------------------------------------------------------------------
# Fit with informative prior
prior_informative <- list(
  rate = list(family = "gamma", hyperparameters = list(shape = 10, rate = 5))
)

fit_informative <- tk_fit(
  data = data,
  censoring_scheme = "complete",
  pdf = pdf_exp,
  cdf = cdf_exp,
  prior_spec = prior_informative,
  initial_values = c(rate = 1),
  loss_function = "sel"
)

# Fit with weakly informative prior
prior_weak <- list(
  rate = list(family = "gamma", hyperparameters = list(shape = 0.1, rate = 0.1))
)

fit_weak <- tk_fit(
  data = data,
  censoring_scheme = "complete",
  pdf = pdf_exp,
  cdf = cdf_exp,
  prior_spec = prior_weak,
  initial_values = c(rate = 1),
  loss_function = "sel"
)

# Compare estimates
data.frame(
  Informative = coef(fit_informative),
  Weak = coef(fit_weak),
  Flat = coef(fit_flat)
)

## -----------------------------------------------------------------------------
sensitivity <- tk_sensitivity(
  fit = fit_informative,
  parameter_name = "rate",
  hyperparameter_name = "shape",
  hyperparameter_values = c(0.1, 0.5, 1, 2, 5, 10)
)

print(sensitivity)
plot(sensitivity)