Tierney-Kadane Approximation: Theory and Derivation

Step 1: Posterior Mode

Define the normalized log-posterior:

\[h(\theta) = \frac{1}{n} [\log L(\theta | x) + \log \pi(\theta)]\]

Find the posterior mode:

\[\hat{\theta} = \arg\max h(\theta)\]

Compute the Hessian at the mode:

\[H(\hat{\theta}) = \frac{\partial^2 h}{\partial \theta \partial \theta^T}\bigg|_{\theta = \hat{\theta}}\]

The approximate posterior covariance is:

\[\Sigma = -H(\hat{\theta})^{-1}\]

Step 2: Modified Function for Expectation

To compute \(E[g(\theta) | x]\), define:

\[h^*(\theta) = h(\theta) + \frac{1}{n} \log g(\theta)\]

Find the mode of the modified function:

\[\theta^* = \arg\max h^*(\theta)\]

Compute the Hessian at this mode:

\[H^*(\theta^*) = \frac{\partial^2 h^*}{\partial \theta \partial \theta^T}\bigg|_{\theta = \theta^*}\]

The modified covariance is:

\[\Sigma^* = -H^*(\theta^*)^{-1}\]

Step 3: TK Approximation Formula

The posterior expectation can be written as:

\[E[g(\theta) | x] = \frac{\int g(\theta) L(\theta | x) \pi(\theta) d\theta}{\int L(\theta | x) \pi(\theta) d\theta} = \frac{\int e^{n h^*(\theta)} d\theta}{\int e^{n h(\theta)} d\theta}\]

Applying the Laplace approximation to both numerator and denominator:

\[E[g(\theta) | x] \approx \frac{e^{n h^*(\theta^*)} \left(\frac{2\pi}{n}\right)^{k/2} |\Sigma^*|^{1/2}}{e^{n h(\hat{\theta})} \left(\frac{2\pi}{n}\right)^{k/2} |\Sigma|^{1/2}}\]

Simplifying:

\[E[g(\theta) | x] \approx \sqrt{\frac{|\Sigma^*|}{|\Sigma|}} \cdot e^{n [h^*(\theta^*) - h(\hat{\theta})]}\]

This is the Tierney-Kadane approximation.

Accuracy

The TK approximation has error of order \(O(n^{-2})\), which is better than the standard Laplace approximation’s \(O(n^{-1})\) error rate.

Requirements

1. Smoothness: \(h(\theta)\) and \(h^*(\theta)\) must be smooth (at least twice differentiable)
2. Unimodality: The posterior should be approximately unimodal
3. Positive definiteness: The Hessian matrices must be positive definite at the modes
4. Sample size: Larger sample sizes improve accuracy

Advantages

1. Computationally efficient: Only requires optimization and Hessian computation
2. No sampling: Unlike MCMC, no random sampling is required
3. Deterministic: Results are reproducible
4. Covariance matrix: Provides approximate posterior covariance directly

Limitations

1. Multimodal posteriors: May fail for highly multimodal posteriors
2. Small samples: Less accurate for very small sample sizes
3. Boundary issues: Can have problems with parameters near boundaries
4. Numerical stability: Requires careful numerical implementation

Function Structure

TKApprox implements the TK approximation through the following functions:

1. tk_posterior(): Constructs \(h(\theta)\) from likelihood and prior
2. tk_mode(): Finds \(\hat{\theta}\) via numerical optimization
3. tk_hessian(): Computes \(H(\hat{\theta})\) and \(\Sigma\)
4. tk_expectation(): Implements the full TK approximation for \(E[g(\theta) | x]\)

Example: Computing Posterior Mean

# Define exponential distribution
pdf_exp <- function(x, param) dexp(x, rate = param)
cdf_exp <- function(x, param) pexp(x, rate = param)

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

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

# Construct log-posterior
log_post <- tk_posterior(data, "complete", pdf_exp, cdf_exp, prior_spec)

# Find posterior mode
mode_result <- tk_mode(log_post, initial_values = c(rate = 1))

# Compute Hessian and covariance
hessian_result <- tk_hessian(log_post, mode_result$mode)

# Compute posterior mean using TK approximation
g_fn <- function(param) param[1]  # g(θ) = θ for posterior mean
tk_result <- tk_expectation(log_post, g_fn, mode_result, hessian_result)

# Compare with direct fit
fit <- tk_fit(data, "complete", pdf_exp, cdf_exp, prior_spec,
              initial_values = c(rate = 1), loss_function = "sel")

data.frame(
  TK_approximation = tk_result$expectation,
  Direct_fit = coef(fit),
  Posterior_mode = mode_result$mode
)

##      TK_approximation Direct_fit Posterior_mode
## rate         1.862275   1.862275       1.777306

Squared Error Loss (Posterior Mean)

\[g(\theta) = \theta\]

\[\hat{\theta}_{SEL} = E[\theta | x]\]

LINEX Loss

\[g(\theta) = \exp(-c\theta)\]

\[\hat{\theta}_{LINEX} = -\frac{1}{c} \log E[\exp(-c\theta) | x]\]

General Entropy Loss

\[g(\theta) = \theta^{-q}\]

\[\hat{\theta}_{GEL} = [E(\theta^{-q} | x)]^{-1/q}\]

Precautionary Loss

\[g(\theta) = \theta^2\]

\[\hat{\theta}_{Precautionary} = \sqrt{E(\theta^2 | x)}\]

Weighted Squared Error Loss

\[g(\theta) = \theta^{-1}\]

\[\hat{\theta}_{WSEL} = \frac{1}{E(\theta^{-1} | x)}\]

Optimization

TKApprox supports multiple optimization methods: - BFGS: Quasi-Newton method (default) - L-BFGS-B: BFGS with box constraints - Nelder-Mead: Simplex method - nlminb: R’s general-purpose optimizer - maxLik: Maximum likelihood package - trust: Trust region methods

Automatic fallback is implemented if the primary method fails.

Numerical Differentiation

The package uses numDeriv for numerical differentiation: - Richardson extrapolation: Default, high accuracy - Simple finite differences: Faster but less accurate

Users can optionally provide analytic derivatives for improved accuracy and speed.

Numerical Stability

Several safeguards are implemented: - Checking for positive definite Hessians - Handling overflow/underflow in likelihood evaluation - Detecting NaN/Inf propagation - Warning on failed optimizations - Alternative optimizer fallback

vs. MCMC

vs. Variational Inference

vs. Importance Sampling