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

## -----------------------------------------------------------------------------
# Define log-normal PDF
pdf_lognormal <- function(x, param) {
  dlnorm(x, meanlog = param[1], sdlog = param[2])
}

# Define log-normal CDF
cdf_lognormal <- function(x, param) {
  plnorm(x, meanlog = param[1], sdlog = param[2])
}

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

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

# Fit the model
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)

## -----------------------------------------------------------------------------
# Define Pareto PDF
pdf_pareto <- function(x, param) {
  xm <- param[1]  # scale parameter (minimum)
  alpha <- param[2]  # shape parameter
  ifelse(x >= xm, (alpha * xm^alpha) / (x^(alpha + 1)), 0)
}

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

# Specify priors
prior_spec <- list(
  xm = list(family = "gamma", hyperparameters = list(shape = 2, rate = 1)),
  alpha = list(family = "gamma", hyperparameters = list(shape = 2, rate = 1))
)

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

# Fit the model
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)

## -----------------------------------------------------------------------------
# Define Burr Type XII PDF
pdf_burr <- function(x, param) {
  c <- param[1]  # shape parameter 1
  k <- param[2]  # shape parameter 2
  lambda <- param[3]  # scale parameter
  (c * k / lambda) * (x / lambda)^(c - 1) / (1 + (x / lambda)^c)^(k + 1)
}

# Define Burr Type XII CDF
cdf_burr <- function(x, param) {
  c <- param[1]
  k <- param[2]
  lambda <- param[3]
  1 - 1 / (1 + (x / lambda)^c)^k
}

# Specify priors
prior_spec <- list(
  c = list(family = "gamma", hyperparameters = list(shape = 2, rate = 1)),
  k = list(family = "gamma", hyperparameters = list(shape = 2, rate = 1)),
  lambda = list(family = "gamma", hyperparameters = list(shape = 2, rate = 1))
)

# Generate Burr data (using approximation)
set.seed(123)
data <- rburr <- function(n, c, k, lambda) {
  u <- runif(n)
  lambda * (u^(-1/k) - 1)^(-1/c)
}
data <- rburr(20, c = 2, k = 1, lambda = 1)

# Fit the model
fit <- tk_fit(
  data = data,
  censoring_scheme = "complete",
  pdf = pdf_burr,
  cdf = cdf_burr,
  prior_spec = prior_spec,
  initial_values = c(c = 1.5, k = 0.8, lambda = 1),
  loss_function = "sel"
)

summary(fit)

## -----------------------------------------------------------------------------
# Define Poisson PMF
pmf_poisson <- function(x, param) {
  dpois(x, lambda = param[1])
}

# For discrete distributions, CDF is still needed for censoring
cdf_poisson <- function(x, param) {
  ppois(x, lambda = param[1])
}

# Specify prior
prior_spec <- list(
  lambda = list(family = "gamma", hyperparameters = list(shape = 2, rate = 1))
)

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

# Fit the model
fit <- tk_fit(
  data = data,
  censoring_scheme = "complete",
  pdf = pmf_poisson,  # Use pmf as pdf for discrete
  cdf = cdf_poisson,
  prior_spec = prior_spec,
  initial_values = c(lambda = 2),
  loss_function = "sel"
)

summary(fit)

## -----------------------------------------------------------------------------
# Define Negative Binomial PMF
pmf_nbinom <- function(x, param) {
  size <- param[1]
  mu <- param[2]
  dnbinom(x, size = size, mu = mu)
}

# Define Negative Binomial CDF
cdf_nbinom <- function(x, param) {
  size <- param[1]
  mu <- param[2]
  pnbinom(x, size = size, mu = mu)
}

# Specify priors
prior_spec <- list(
  size = list(family = "gamma", hyperparameters = list(shape = 2, rate = 1)),
  mu = list(family = "gamma", hyperparameters = list(shape = 2, rate = 1))
)

# Generate Negative Binomial data
set.seed(123)
data <- rnbinom(20, size = 5, mu = 3)

# Fit the model
fit <- tk_fit(
  data = data,
  censoring_scheme = "complete",
  pdf = pmf_nbinom,
  cdf = cdf_nbinom,
  prior_spec = prior_spec,
  initial_values = c(size = 4, mu = 2.5),
  loss_function = "sel"
)

summary(fit)

## -----------------------------------------------------------------------------
# Custom log-likelihood for a mixture model
loglik_mixture <- function(param, data) {
  p <- param[1]  # mixing proportion
  lambda1 <- param[2]  # rate for component 1
  lambda2 <- param[3]  # rate for component 2
  
  # Enforce constraints
  if (p <= 0 || p >= 1 || lambda1 <= 0 || lambda2 <= 0) {
    return(-Inf)
  }
  
  # Log-likelihood
  ll <- sum(log(p * dexp(data, rate = lambda1) + (1 - p) * dexp(data, rate = lambda2)))
  
  if (!is.finite(ll)) {
    return(-Inf)
  }
  
  ll
}

# Specify priors
prior_spec <- list(
  p = list(family = "beta", hyperparameters = list(shape1 = 2, shape2 = 2)),
  lambda1 = list(family = "gamma", hyperparameters = list(shape = 2, rate = 1)),
  lambda2 = list(family = "gamma", hyperparameters = list(shape = 2, rate = 1))
)

# Generate mixture data
set.seed(123)
data <- c(rexp(10, rate = 1), rexp(10, rate = 5))

# Fit the model
fit <- tk_fit(
  data = data,
  censoring_scheme = "complete",
  pdf = NULL,  # Not needed when providing loglik
  cdf = NULL,
  prior_spec = prior_spec,
  loglik = loglik_mixture,
  initial_values = c(p = 0.5, lambda1 = 1, lambda2 = 5),
  loss_function = "sel"
)

summary(fit)