Defining Custom Distributions

The hardware and bandwidth for this mirror is donated by METANET, the Webhosting and Full Service-Cloud Provider.
If you wish to report a bug, or if you are interested in having us mirror your free-software or open-source project, please feel free to contact us at mirror[@]metanet.ch.

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:

PDF/PMF function: pdf(x, param) or pmf(x, param) for discrete distributions
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

Parameter ordering: Be consistent with parameter ordering in PDF and CDF functions
Boundary handling: Ensure your functions handle boundary cases gracefully
Numerical stability: Use log-scale computations when possible to avoid overflow/underflow
Parameter constraints: Use bounds in tk_fit() to keep parameters in valid ranges
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

These binaries (installable software) and packages are in development.
They may not be fully stable and should be used with caution. We make no claims about them.