---
title: "Defining Custom Distributions"
author: "Your Name"
date: "`r Sys.Date()`"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Defining Custom Distributions}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
library(TKApprox)
```

## Introduction

TKApprox is designed to work with any user-specified univariate probability distribution. This vignette shows how to define custom distributions for use with the package.

## Required Functions

To use a distribution with TKApprox, you typically need to provide:

1. **PDF/PMF function**: `pdf(x, param)` or `pmf(x, param)` for discrete distributions
2. **CDF function**: `cdf(x, param)` (required for censored data)

The parameter `param` is always a numeric vector containing all distribution parameters.

## Continuous Distributions

### Example 1: Log-Normal Distribution

```{r}
# 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)
```

### Example 2: Pareto Distribution

```{r}
# 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)
```

### Example 3: Burr Type XII Distribution

```{r}
# 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)
```

## Discrete Distributions

### Example 4: Poisson Distribution

```{r}
# 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)
```

### Example 5: Negative Binomial Distribution

```{r}
# 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)
```

## Using Custom Log-Likelihood Functions

For complex models, you can provide a custom log-likelihood function directly instead of relying on the automatic construction from PDF/CDF.

```{r}
# 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)
```

## Tips for Defining Distributions

1. **Parameter ordering**: Be consistent with parameter ordering in PDF and CDF functions
2. **Boundary handling**: Ensure your functions handle boundary cases gracefully
3. **Numerical stability**: Use log-scale computations when possible to avoid overflow/underflow
4. **Parameter constraints**: Use bounds in `tk_fit()` to keep parameters in valid ranges
5. **Initial values**: Good initial values are crucial for convergence in complex models

## Next Steps

- See "Censoring Schemes" for how to handle censored data with custom distributions
- See "Prior Specification" for advanced prior modeling techniques