Introduction to TKApprox

Introduction

TKApprox provides a distribution-independent framework for Bayesian estimation of arbitrary univariate probability models using the Tierney-Kadane approximation. This vignette provides a gentle introduction to the package’s main features and workflow.

Basic Workflow

The typical workflow in TKApprox follows these steps:

Define your probability distribution (PDF/PMF and CDF)
Specify prior distributions for parameters
Provide your data and censoring scheme
Call tk_fit() to perform Bayesian estimation
Examine results using S3 methods and visualization

A Simple Example: Exponential Distribution

Let’s start with a simple example using the exponential distribution.

Step 1: Define the Distribution

# Define the exponential PDF
pdf_exp <- function(x, param) {
  dexp(x, rate = param)
}

# Define the exponential CDF
cdf_exp <- function(x, param) {
  pexp(x, rate = param)
}

Step 2: Specify Prior

# Gamma prior for the rate parameter
prior_spec <- list(
  rate = list(family = "gamma", hyperparameters = list(shape = 2, rate = 1))
)

Step 3: Generate Data

set.seed(123)
data <- rexp(20, rate = 1.5)

Step 4: Fit the Model

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"
)

Step 5: Examine Results

# Print basic results
print(fit)

## Tierney-Kadane Bayesian Estimation
## ===================================
## 
## Censoring scheme: complete 
## Sample size: 20 
## Loss function: sel 
## Optimization method: nlminb 
## Convergence: 0 
## Iterations: 7 
## Execution time: 0.0377 seconds
## 
## Posterior mode (MAP):
##     rate 
## 1.777306 
## 
## Bayes estimates (sel):
##     rate 
## 1.862275 
## 
## Standard errors:
##    rate 
## 0.38784 
## 
## 95% Credible intervals (normal approximation):
##         2.5 %   97.5 %
## rate 1.102123 2.622428
## 
## Log-posterior at mode: -0.4461 
## Log-likelihood at mode: -7.7207

# Detailed summary
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: nlminb 
## Convergence code: 0 
## Iterations: 7 
## Gradient norm: 0 
## Execution time: 0.0377 seconds
## 
## Parameter Estimates:
## --------------------
##  Parameter Posterior_Mode Bayes_Estimate Std_Error CI_Lower CI_Upper
##       rate       1.777306       1.862275   0.38784 1.102123 2.622428
## 
## Model Fit Statistics:
## ---------------------
## Log-posterior at mode: -0.4461 
## Log-likelihood at mode: -7.7207 
## Prior contribution: -1.2022 
## 
## Posterior Covariance Matrix:
## ---------------------------
##         rate
## rate 0.15042

# Extract Bayes estimates
coef(fit)

##     rate 
## 1.862275

# Extract covariance matrix
vcov(fit)

##           rate
## rate 0.1504198

# Model comparison statistics
print_model_comparison(fit)

## 
## === Model Comparison Statistics ===
## 
##                Statistic      Value
##           Log-Likelihood -7.7207187
##  Negative Log-Likelihood  7.7207187
##                      AIC 17.4414374
##                      BIC 18.4371697
##                     CAIC 19.4371697
##                     HQIC 17.6358148
##                      DIC 17.4414374
##   Expected Log-Posterior -0.4461463
##     Number of Parameters  1.0000000
##              Sample Size 20.0000000

Visualization

# Plot all diagnostics
plot(fit)

# Or select specific plots
plot(fit, which = 1)  # Posterior approximation

plot(fit, which = 2)  # Likelihood surface

plot(fit, which = 3)  # Prior vs posterior

Different Loss Functions

TKApprox supports multiple Bayesian loss functions:

Squared Error Loss (Posterior Mean)

fit_sel <- 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"
)
coef(fit_sel)

##     rate 
## 1.862275

LINEX Loss

fit_linex <- tk_fit(
  data = data,
  censoring_scheme = "complete",
  pdf = pdf_exp,
  cdf = cdf_exp,
  prior_spec = prior_spec,
  initial_values = c(rate = 1),
  loss_function = "linex",
  loss_params = list(c = 0.5)
)
coef(fit_linex)

##     rate 
## 1.823621

General Entropy Loss

fit_gel <- tk_fit(
  data = data,
  censoring_scheme = "complete",
  pdf = pdf_exp,
  cdf = cdf_exp,
  prior_spec = prior_spec,
  initial_values = c(rate = 1),
  loss_function = "gel",
  loss_params = list(q = 0.5)
)
coef(fit_gel)

##     rate 
## 1.798937

Censored Data

TKApprox handles various censoring schemes. Here’s an example with right-censored data:

# Create right-censored data
status <- c(1, 1, 0, 1, 0, 1, 1, 0, 1, 1)  # 1 = observed, 0 = censored

fit_censored <- tk_fit(
  data = data,
  censoring_scheme = "right-censored",
  pdf = pdf_exp,
  cdf = cdf_exp,
  prior_spec = prior_spec,
  initial_values = c(rate = 1),
  loss_function = "sel",
  status = status
)

summary(fit_censored)

## 
## === Tierney-Kadane Bayesian Estimation Summary ===
## 
## Model Information:
## -----------------
## Censoring scheme: right-censored 
## Sample size: 20 
## Number of parameters: 1 
## Loss function: sel 
## 
## Optimization Results:
## --------------------
## Method: BFGS 
## Convergence code: 0 
## Iterations: 1 
## Gradient norm: 0 
## Execution time: 0.01 seconds
## 
## Parameter Estimates:
## --------------------
##  Parameter Posterior_Mode Bayes_Estimate Std_Error  CI_Lower CI_Upper
##       rate              1              1 0.2236068 0.5617387 1.438261
## 
## Model Fit Statistics:
## ---------------------
## Log-posterior at mode: -1e+10 
## Log-likelihood at mode: -Inf 
## Prior contribution: -1 
## 
## Posterior Covariance Matrix:
## ---------------------------
##      rate
## rate 0.05

Prior Sensitivity Analysis

Examine how sensitive your estimates are to prior hyperparameters:

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

print(sensitivity)

## Prior Sensitivity Analysis
## ==========================
## Parameter: rate 
## Hyperparameter: shape 
## Loss function: sel 
## Number of hyperparameter values tested: 5 
## 
## Results:
##  hyperparameter_value log_posterior log_likelihood convergence iterations
##                   0.5             0              0           0          0
##                   1.0             0              0           0          0
##                   2.0             0              0           0          0
##                   5.0             0              0           0          0
##                  10.0             0              0           0          0
##  estimate_rate se_rate
##             NA      NA
##             NA      NA
##             NA      NA
##             NA      NA
##             NA      NA

plot(sensitivity)

## True parameter not available; cannot compute risk.

Multi-Parameter Models

TKApprox works with models of any dimensionality. Here’s a two-parameter example with the Weibull distribution:

# Define 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
prior_spec_weibull <- list(
  shape = list(family = "gamma", hyperparameters = list(shape = 2, rate = 1)),
  scale = list(family = "gamma", hyperparameters = list(shape = 2, rate = 1))
)

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

# Fit the model
fit_weibull <- tk_fit(
  data = data_weibull,
  censoring_scheme = "complete",
  pdf = pdf_weibull,
  cdf = cdf_weibull,
  prior_spec = prior_spec_weibull,
  initial_values = c(shape = 1.5, scale = 1),
  loss_function = "sel"
)

summary(fit_weibull)

## 
## === 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: 1e-06 
## Execution time: 0.1372 seconds
## 
## Parameter Estimates:
## --------------------
##  Parameter Posterior_Mode Bayes_Estimate Std_Error CI_Lower CI_Upper
##      shape      1.7566163      1.7595973 0.3066255 1.158622 2.360572
##      scale      0.9124234      0.9476932 0.1210758 0.710389 1.184997
## 
## Model Fit Statistics:
## ---------------------
## Log-posterior at mode: -0.6901 
## Log-likelihood at mode: -11.6046 
## Prior contribution: -2.1973 
## 
## Posterior Covariance Matrix:
## ---------------------------
##          shape    scale
## shape 0.094019 0.011196
## scale 0.011196 0.014659

plot(fit_weibull, which = 2)  # Likelihood surface for 2-parameter model

Next Steps

See “Defining Custom Distributions” for more complex distribution examples
See “Censoring Schemes” for detailed coverage of all censoring types
See “Loss Functions” for more information on Bayesian loss functions
See “Prior Specification” for advanced prior modeling