Defining Custom Distributions

Your Name

2026-07-31

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

# 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)
## 
## === Tierney-Kadane Bayesian Estimation Summary ===
## 
## Model Information:
## -----------------
## Censoring scheme: complete 
## Sample size: 20 
## Number of parameters: 2 
## Loss function: sel 
## 
## Optimization Results:
## --------------------
## Method: nlminb 
## Convergence code: 0 
## Iterations: 6 
## Gradient norm: 2e-06 
## Execution time: 0.1695 seconds
## 
## Parameter Estimates:
## --------------------
##  Parameter Posterior_Mode Bayes_Estimate Std_Error   CI_Lower  CI_Upper
##    meanlog     0.07000461      0.0960029 0.1067846 -0.1132912 0.3052970
##      sdlog     0.48030034      0.5255941 0.0764787  0.3756986 0.6754896
## 
## Model Fit Statistics:
## ---------------------
## Log-posterior at mode: -0.8502 
## Log-likelihood at mode: -14.8683 
## Prior contribution: -2.135 
## 
## Posterior Covariance Matrix:
## ---------------------------
##           meanlog     sdlog
## meanlog  0.011403 -0.000019
## sdlog   -0.000019  0.005849

Example 2: Pareto Distribution

# 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)
## 
## === Tierney-Kadane Bayesian Estimation Summary ===
## 
## Model Information:
## -----------------
## Censoring scheme: complete 
## Sample size: 20 
## Number of parameters: 2 
## Loss function: sel 
## 
## Optimization Results:
## --------------------
## Method: BFGS 
## Convergence code: 0 
## Iterations: 93 
## Gradient norm: 1.814221 
## Execution time: 0.4392 seconds
## 
## Parameter Estimates:
## --------------------
##  Parameter Posterior_Mode Bayes_Estimate Std_Error CI_Lower CI_Upper
##         xm       1.021130       1.021130 0.2236068 0.582869 1.459392
##      alpha       1.852707       1.852707 0.2236068 1.414446 2.290968
## 
## Model Fit Statistics:
## ---------------------
## Log-posterior at mode: -1.152 
## Log-likelihood at mode: -20.8028 
## Prior contribution: -2.2363 
## 
## Posterior Covariance Matrix:
## ---------------------------
##         xm alpha
## xm    0.05  0.00
## alpha 0.00  0.05

Example 3: Burr Type XII Distribution

# 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)
## 
## === Tierney-Kadane Bayesian Estimation Summary ===
## 
## Model Information:
## -----------------
## Censoring scheme: complete 
## Sample size: 20 
## Number of parameters: 3 
## Loss function: sel 
## 
## Optimization Results:
## --------------------
## Method: nlminb 
## Convergence code: 0 
## Iterations: 13 
## Gradient norm: 0 
## Execution time: 0.5093 seconds
## 
## Parameter Estimates:
## --------------------
##  Parameter Posterior_Mode Bayes_Estimate Std_Error  CI_Lower CI_Upper
##          c       1.705931       1.605225 0.4002586 0.8207327 2.389717
##          k       1.216901       2.007805 0.6891115 0.6571714 3.358439
##     lambda       1.375214       2.258235 0.7132158 0.8603579 3.656113
## 
## Model Fit Statistics:
## ---------------------
## Log-posterior at mode: -1.67 
## Log-likelihood at mode: -30.1503 
## Prior contribution: -3.249 
## 
## Posterior Covariance Matrix:
## ---------------------------
##                c         k    lambda
## c       0.160207 -0.178914 -0.173239
## k      -0.178914  0.474875  0.449851
## lambda -0.173239  0.449851  0.508677

Discrete Distributions

Example 4: Poisson Distribution

# 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)
## 
## === Tierney-Kadane Bayesian Estimation Summary ===
## 
## Model Information:
## -----------------
## Censoring scheme: complete 
## Sample size: 20 
## Number of parameters: 1 
## Loss function: sel 
## 
## Optimization Results:
## --------------------
## Method: L-BFGS-B 
## Convergence code: 0 
## Iterations: 6 
## Gradient norm: 0 
## Execution time: 0.0384 seconds
## 
## Parameter Estimates:
## --------------------
##  Parameter Posterior_Mode Bayes_Estimate Std_Error CI_Lower CI_Upper
##     lambda       3.142857       3.190536  0.386859 2.432307 3.948766
## 
## Model Fit Statistics:
## ---------------------
## Log-posterior at mode: -2.1304 
## Log-likelihood at mode: -40.6113 
## Prior contribution: -1.9977 
## 
## Posterior Covariance Matrix:
## ---------------------------
##         lambda
## lambda 0.14966

Example 5: Negative Binomial Distribution

# 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)
## 
## === Tierney-Kadane Bayesian Estimation Summary ===
## 
## Model Information:
## -----------------
## Censoring scheme: complete 
## Sample size: 20 
## Number of parameters: 2 
## Loss function: sel 
## 
## Optimization Results:
## --------------------
## Method: nlminb 
## Convergence code: 0 
## Iterations: 13 
## Gradient norm: 0 
## Execution time: 0.1906 seconds
## 
## Parameter Estimates:
## --------------------
##  Parameter Posterior_Mode Bayes_Estimate Std_Error CI_Lower CI_Upper
##       size       3.376409       4.089861  1.399050 1.347774 6.831949
##         mu       3.187384       3.316787  0.525078 2.287653 4.345921
## 
## Model Fit Statistics:
## ---------------------
## Log-posterior at mode: -2.3119 
## Log-likelihood at mode: -42.0505 
## Prior contribution: -4.1878 
## 
## Posterior Covariance Matrix:
## ---------------------------
##          size       mu
## size 1.957341 0.052986
## mu   0.052986 0.275707

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.

# 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)
## 
## === Tierney-Kadane Bayesian Estimation Summary ===
## 
## Model Information:
## -----------------
## Censoring scheme: complete 
## Sample size: 20 
## Number of parameters: 3 
## Loss function: sel 
## 
## Optimization Results:
## --------------------
## Method: nlminb 
## Convergence code: 0 
## Iterations: 15 
## Gradient norm: 0 
## Execution time: 0.3327 seconds
## 
## Parameter Estimates:
## --------------------
##  Parameter Posterior_Mode Bayes_Estimate Std_Error  CI_Lower  CI_Upper
##          p      0.2811991       0.445935 0.1857383 0.0818947 0.8099753
##    lambda1      1.0067579       1.534428 0.5418351 0.4724511 2.5964058
##    lambda2      3.7701084       4.329415 1.2741434 1.8321400 6.8266902
## 
## Model Fit Statistics:
## ---------------------
## Log-posterior at mode: -0.1907 
## Log-likelihood at mode: -0.5645 
## Prior contribution: -3.2501 
## 
## Posterior Covariance Matrix:
## ---------------------------
##                p  lambda1  lambda2
## p       0.034499 0.053345 0.071737
## lambda1 0.053345 0.293585 0.056169
## lambda2 0.071737 0.056169 1.623441

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