Prior Specification

Introduction

TKApprox provides flexible prior specification options, from standard conjugate priors to fully custom prior functions. This vignette covers all available prior families and how to use them effectively.

Standard Prior Families

Gamma Prior

The Gamma prior is commonly used for positive parameters like rates and scales.

Parameters: shape (α), rate (β)

PDF: \(f(x) = \frac{\beta^\alpha}{\Gamma(\alpha)} x^{\alpha-1} e^{-\beta x}\)

# Gamma prior for exponential rate parameter
pdf_exp <- function(x, param) dexp(x, rate = param)
cdf_exp <- function(x, param) pexp(x, rate = param)

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

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

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

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.033 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

Normal Prior

The Normal prior is used for parameters that can take any real value.

Parameters: mean (μ), standard deviation (σ)

PDF: \(f(x) = \frac{1}{\sqrt{2\pi}\sigma} \exp\left(-\frac{(x-\mu)^2}{2\sigma^2}\right)\)

# Normal prior for log-normal meanlog parameter
pdf_lognormal <- function(x, param) dlnorm(x, meanlog = param[1], sdlog = param[2])
cdf_lognormal <- function(x, param) plnorm(x, meanlog = param[1], sdlog = param[2])

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

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

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.1393 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

Beta Prior

The Beta prior is used for parameters bounded between 0 and 1.

Parameters: shape1 (α), shape2 (β)

PDF: \(f(x) = \frac{x^{\alpha-1}(1-x)^{\beta-1}}{B(\alpha,\beta)}\)

# Beta prior for probability parameter
pdf_bernoulli <- function(x, param) {
  p <- param[1]
  ifelse(x == 1, p, 1 - p)
}

cdf_bernoulli <- function(x, param) {
  p <- param[1]
  ifelse(x == 0, 1 - p, 1)
}

prior_spec <- list(
  p = list(family = "beta", hyperparameters = list(shape1 = 2, shape2 = 2))
)

# Bernoulli data
set.seed(123)
data <- rbinom(20, size = 1, prob = 0.6)

fit <- tk_fit(
  data = data,
  censoring_scheme = "complete",
  pdf = pdf_bernoulli,
  cdf = cdf_bernoulli,
  prior_spec = prior_spec,
  initial_values = c(p = 0.5),
  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: nlminb 
## Convergence code: 0 
## Iterations: 5 
## Gradient norm: 0 
## Execution time: 0.0621 seconds
## 
## Parameter Estimates:
## --------------------
##  Parameter Posterior_Mode Bayes_Estimate Std_Error  CI_Lower  CI_Upper
##          p      0.5909091      0.5824012 0.1048236 0.3769508 0.7878516
## 
## Model Fit Statistics:
## ---------------------
## Log-posterior at mode: -0.6546 
## Log-likelihood at mode: -13.4637 
## Prior contribution: 0.3718 
## 
## Posterior Covariance Matrix:
## ---------------------------
##          p
## p 0.010988

Uniform Prior

The Uniform prior represents a non-informative prior over a bounded interval.

Parameters: lower (a), upper (b)

PDF: \(f(x) = \frac{1}{b-a}\) for \(a \leq x \leq b\)

# Uniform prior for Weibull shape parameter
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])

prior_spec <- list(
  shape = list(family = "uniform", hyperparameters = list(lower = 0.1, upper = 10)),
  scale = list(family = "gamma", hyperparameters = list(shape = 2, rate = 1))
)

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

fit <- tk_fit(
  data = data,
  censoring_scheme = "complete",
  pdf = pdf_weibull,
  cdf = cdf_weibull,
  prior_spec = prior_spec,
  initial_values = c(shape = 1.5, scale = 1),
  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: 7 
## Gradient norm: 0 
## Execution time: 0.1668 seconds
## 
## Parameter Estimates:
## --------------------
##  Parameter Posterior_Mode Bayes_Estimate Std_Error CI_Lower CI_Upper
##      shape       1.799199      1.8021470 0.3174707 1.179916 2.424378
##      scale       0.917478      0.9514555 0.1194147 0.717407 1.185504
## 
## Model Fit Statistics:
## ---------------------
## Log-posterior at mode: -0.7446 
## Log-likelihood at mode: -11.5959 
## Prior contribution: -3.2961 
## 
## Posterior Covariance Matrix:
## ---------------------------
##          shape    scale
## shape 0.100788 0.011924
## scale 0.011924 0.014260

Exponential Prior

The Exponential prior is a special case of Gamma with shape = 1.

Parameters: rate (λ)

PDF: \(f(x) = \lambda e^{-\lambda x}\)

# Exponential prior for Poisson rate
pdf_poisson <- function(x, param) dpois(x, lambda = param[1])
cdf_poisson <- function(x, param) ppois(x, lambda = param[1])

prior_spec <- list(
  lambda = list(family = "exponential", hyperparameters = list(rate = 1))
)

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

fit <- tk_fit(
  data = data,
  censoring_scheme = "complete",
  pdf = pdf_poisson,
  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.0342 seconds
## 
## Parameter Estimates:
## --------------------
##  Parameter Posterior_Mode Bayes_Estimate Std_Error CI_Lower CI_Upper
##     lambda       3.095238       3.142918  0.383917 2.390454 3.895382
## 
## Model Fit Statistics:
## ---------------------
## Log-posterior at mode: -2.1873 
## Log-likelihood at mode: -40.6513 
## Prior contribution: -3.0952 
## 
## Posterior Covariance Matrix:
## ---------------------------
##          lambda
## lambda 0.147392

Log-Normal Prior

The Log-Normal prior is useful for parameters that are log-normally distributed.

Parameters: meanlog (μ), sdlog (σ)

PDF: \(f(x) = \frac{1}{x\sigma\sqrt{2\pi}} \exp\left(-\frac{(\log x - \mu)^2}{2\sigma^2}\right)\)

# Log-Normal prior for Pareto scale parameter
pdf_pareto <- function(x, param) {
  xm <- param[1]
  alpha <- param[2]
  ifelse(x >= xm, (alpha * xm^alpha) / (x^(alpha + 1)), 0)
}

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

prior_spec <- list(
  xm = list(family = "lognormal", hyperparameters = list(meanlog = 0, sdlog = 0.5)),
  alpha = list(family = "gamma", hyperparameters = list(shape = 2, rate = 1))
)

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

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: 108 
## Gradient norm: 1.717385 
## Execution time: 0.6182 seconds
## 
## Parameter Estimates:
## --------------------
##  Parameter Posterior_Mode Bayes_Estimate Std_Error  CI_Lower CI_Upper
##         xm       1.021213       1.021213 0.2236068 0.5829514 1.459474
##      alpha       1.807477       1.807477 0.2236068 1.3692161 2.245739
## 
## Model Fit Statistics:
## ---------------------
## Log-posterior at mode: -1.1119 
## Log-likelihood at mode: -20.7755 
## Prior contribution: -1.4632 
## 
## Posterior Covariance Matrix:
## ---------------------------
##         xm alpha
## xm    0.05  0.00
## alpha 0.00  0.05

Weibull Prior

The Weibull prior is useful for reliability and survival analysis parameters.

Parameters: shape (k), scale (λ)

PDF: \(f(x) = \frac{k}{\lambda}\left(\frac{x}{\lambda}\right)^{k-1} e^{-(x/\lambda)^k}\)

# Weibull prior for gamma shape parameter
pdf_gamma <- function(x, param) dgamma(x, shape = param[1], rate = param[2])
cdf_gamma <- function(x, param) pgamma(x, shape = param[1], rate = param[2])

prior_spec <- list(
  shape = list(family = "weibull", hyperparameters = list(shape = 2, scale = 1)),
  rate = list(family = "gamma", hyperparameters = list(shape = 2, rate = 1))
)

set.seed(123)
data <- rgamma(20, shape = 2, rate = 1.5)

fit <- tk_fit(
  data = data,
  censoring_scheme = "complete",
  pdf = pdf_gamma,
  cdf = cdf_gamma,
  prior_spec = prior_spec,
  initial_values = c(shape = 1.5, rate = 1),
  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: 0 
## Execution time: 0.1731 seconds
## 
## Parameter Estimates:
## --------------------
##  Parameter Posterior_Mode Bayes_Estimate Std_Error  CI_Lower CI_Upper
##      shape       1.297166       1.396159 0.3085844 0.7913446 2.000973
##       rate       1.119369       1.243792 0.3350343 0.5871373 1.900448
## 
## Model Fit Statistics:
## ---------------------
## Log-posterior at mode: -1.1817 
## Log-likelihood at mode: -21.8977 
## Prior contribution: -1.7359 
## 
## Posterior Covariance Matrix:
## ---------------------------
##          shape     rate
## shape 0.095224 0.079123
## rate  0.079123 0.112248

Inverse Gamma Prior

The Inverse Gamma prior is commonly used for variance parameters.

Parameters: shape (α), scale (β)

PDF: \(f(x) = \frac{\beta^\alpha}{\Gamma(\alpha)} x^{-\alpha-1} e^{-\beta/x}\)

# Inverse Gamma prior for normal variance
pdf_normal <- function(x, param) dnorm(x, mean = param[1], sd = sqrt(param[2]))
cdf_normal <- function(x, param) pnorm(x, mean = param[1], sd = sqrt(param[2]))

prior_spec <- list(
  mean = list(family = "normal", hyperparameters = list(mean = 0, sd = 10)),
  variance = list(family = "invgamma", hyperparameters = list(shape = 2, scale = 1))
)

set.seed(123)
data <- rnorm(20, mean = 0, sd = 2)

fit <- tk_fit(
  data = data,
  censoring_scheme = "complete",
  pdf = pdf_normal,
  cdf = cdf_normal,
  prior_spec = prior_spec,
  initial_values = c(mean = 0, variance = 4),
  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: 11 
## Gradient norm: 0 
## Execution time: 0.1659 seconds
## 
## Parameter Estimates:
## --------------------
##  Parameter Posterior_Mode Bayes_Estimate Std_Error   CI_Lower CI_Upper
##       mean      0.2828456      0.3666614 0.3767192 -0.3716948 1.105018
##   variance      2.8423814      3.4700076 0.7883348  1.9248998 5.015115
## 
## Model Fit Statistics:
## ---------------------
## Log-posterior at mode: -2.4091 
## Log-likelihood at mode: -41.4734 
## Prior contribution: -6.7077 
## 
## Posterior Covariance Matrix:
## ---------------------------
##               mean  variance
## mean      0.141917 -0.000088
## variance -0.000088  0.621472

Independent Priors for Multiple Parameters

For multi-parameter models, you can specify independent priors for each parameter:

# Two-parameter 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 for each parameter
prior_spec <- list(
  shape = list(family = "gamma", hyperparameters = list(shape = 2, rate = 1)),
  scale = list(family = "gamma", hyperparameters = list(shape = 2, rate = 1))
)

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

fit <- tk_fit(
  data = data,
  censoring_scheme = "complete",
  pdf = pdf_weibull,
  cdf = cdf_weibull,
  prior_spec = prior_spec,
  initial_values = c(shape = 1.5, scale = 1),
  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: 1e-06 
## Execution time: 0.143 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

Custom Prior Functions

You can also specify a custom prior function directly:

# Custom prior function
custom_logprior <- function(param) {
  # Example: hierarchical prior
  # param[1] = theta, param[2] = hyperparameter
  theta <- param[1]
  hyper <- param[2]
  
  # Prior for theta given hyper
  log_prior_theta <- dnorm(theta, mean = 0, sd = hyper, log = TRUE)
  
  # Prior for hyper
  log_prior_hyper <- dgamma(hyper, shape = 2, rate = 1, log = TRUE)
  
  log_prior_theta + log_prior_hyper
}

# Use custom prior in tk_fit
fit <- tk_fit(
  data = data,
  censoring_scheme = "complete",
  pdf = pdf_exp,
  cdf = cdf_exp,
  prior_spec = custom_logprior,
  initial_values = c(rate = 1),
  loss_function = "sel"
)

Non-Informative (Flat) Priors

To use a non-informative flat prior, simply set prior_spec = NULL:

fit_flat <- tk_fit(
  data = data,
  censoring_scheme = "complete",
  pdf = pdf_exp,
  cdf = cdf_exp,
  prior_spec = NULL,  # Flat prior
  initial_values = c(rate = 1),
  loss_function = "sel"
)

summary(fit_flat)

## 
## === 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: 5 
## Gradient norm: 0 
## Execution time: 0.0282 seconds
## 
## Parameter Estimates:
## --------------------
##  Parameter Posterior_Mode Bayes_Estimate Std_Error CI_Lower CI_Upper
##       rate       1.231553       1.293387 0.2753835 0.753645 1.833129
## 
## Model Fit Statistics:
## ---------------------
## Log-posterior at mode: -0.7917 
## Log-likelihood at mode: -15.8345 
## Prior contribution: 0 
## 
## Posterior Covariance Matrix:
## ---------------------------
##          rate
## rate 0.075836

Prior Sensitivity Analysis

It’s important to check how sensitive your results are to prior specifications:

# Fit with informative prior
prior_informative <- list(
  rate = list(family = "gamma", hyperparameters = list(shape = 10, rate = 5))
)

fit_informative <- tk_fit(
  data = data,
  censoring_scheme = "complete",
  pdf = pdf_exp,
  cdf = cdf_exp,
  prior_spec = prior_informative,
  initial_values = c(rate = 1),
  loss_function = "sel"
)

# Fit with weakly informative prior
prior_weak <- list(
  rate = list(family = "gamma", hyperparameters = list(shape = 0.1, rate = 0.1))
)

fit_weak <- tk_fit(
  data = data,
  censoring_scheme = "complete",
  pdf = pdf_exp,
  cdf = cdf_exp,
  prior_spec = prior_weak,
  initial_values = c(rate = 1),
  loss_function = "sel"
)

# Compare estimates
data.frame(
  Informative = coef(fit_informative),
  Weak = coef(fit_weak),
  Flat = coef(fit_flat)
)

##      Informative     Weak     Flat
## rate    1.412587 1.230403 1.293387

Systematic Prior Sensitivity

Use tk_sensitivity() for systematic examination of prior hyperparameters:

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

print(sensitivity)

## Prior Sensitivity Analysis
## ==========================
## Parameter: rate 
## Hyperparameter: shape 
## Loss function: sel 
## Number of hyperparameter values tested: 6 
## 
## Results:
##  hyperparameter_value log_posterior log_likelihood convergence iterations
##                   0.1             0              0           0          0
##                   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
##             NA      NA

plot(sensitivity)

## True parameter not available; cannot compute risk.

Choosing Prior Hyperparameters

Conjugate Priors

For common distributions, conjugate priors provide computational advantages:

Exponential likelihood + Gamma prior: Conjugate
Normal likelihood (known variance) + Normal prior: Conjugate
Normal likelihood (known mean) + Inverse Gamma prior: Conjugate
Binomial likelihood + Beta prior: Conjugate
Poisson likelihood + Gamma prior: Conjugate

Weakly Informative Priors

When you have little prior information, use weakly informative priors:

Gamma(0.1, 0.1): Very weak prior for positive parameters
Normal(0, 100): Very weak prior for real-valued parameters
Beta(1, 1): Uniform prior for probabilities

Informative Priors

When you have strong prior information (e.g., from previous studies):

Use hyperparameters that reflect your prior knowledge
Consider the prior variance relative to expected data information
Always perform sensitivity analysis

Tips for Prior Specification

Check prior support: Ensure the prior support matches the parameter domain
Scale matters: Consider the scale of your parameters when setting hyperparameters
Sensitivity analysis: Always check how sensitive results are to prior choices
Conjugate when possible: Use conjugate priors for computational efficiency
Document choices: Document your prior specification rationale

Next Steps

See “Loss Functions” for information on Bayesian estimation methods
See “Simulation Studies” for comparing prior specifications in simulation