---
title: "Tierney-Kadane Approximation: Theory and Derivation"
author: "Your Name"
date: "`r Sys.Date()`"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Tierney-Kadane Approximation: Theory and Derivation}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
library(TKApprox)
```

## Introduction

The Tierney-Kadane (TK) approximation is a powerful method for approximating posterior expectations in Bayesian analysis. This vignette provides the theoretical foundation and derivation of the TK approximation as implemented in TKApprox.

## Problem Statement

Given data $x$ and parameter vector $\theta$, we want to compute posterior expectations of the form:

$$E[g(\theta) | x] = \int g(\theta) p(\theta | x) d\theta$$

where $p(\theta | x)$ is the posterior distribution:

$$p(\theta | x) = \frac{L(\theta | x) \pi(\theta)}{\int L(\theta | x) \pi(\theta) d\theta}$$

Here, $L(\theta | x)$ is the likelihood and $\pi(\theta)$ is the prior.

For many models, these integrals are intractable and require numerical methods. The Tierney-Kadane approximation provides an efficient analytical approximation.

## Laplace Approximation Background

The TK approximation builds on the Laplace approximation, which approximates integrals of the form:

$$I = \int e^{n h(\theta)} d\theta$$

where $n$ is the sample size and $h(\theta)$ is a smooth function.

The Laplace approximation states:

$$I \approx e^{n h(\hat{\theta})} \left(\frac{2\pi}{n}\right)^{k/2} |\Sigma|^{1/2}$$

where:
- $\hat{\theta} = \arg\max h(\theta)$ is the mode of $h$
- $\Sigma = -[H(\hat{\theta})]^{-1}$ is the negative inverse Hessian at the mode
- $k$ is the dimension of $\theta$

## Tierney-Kadane 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.

## Properties of the TK 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

## Implementation in TKApprox

### 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

```{r}
# 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
)
```

## Connection to Loss Functions

The TK approximation framework naturally extends to various Bayesian loss functions through appropriate choice of $g(\theta)$:

### 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)}$$

## Numerical Considerations

### 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

## Comparison with Other Methods

### vs. MCMC

| Aspect | TK Approximation | MCMC |
|--------|------------------|------|
| Speed | Fast | Slow |
| Deterministic | Yes | No |
| Multimodal | Limited | Handles well |
| Small samples | Less accurate | Accurate |
| Implementation | Simple | Complex |
| Reproducibility | Perfect | Requires seed |

### vs. Variational Inference

| Aspect | TK Approximation | VI |
|--------|------------------|-----|
| Accuracy | High (large n) | Variable |
| Speed | Fast | Fast |
| Theoretical guarantees | Asymptotic | Optimization-based |
| Implementation | Simple | Complex |

### vs. Importance Sampling

| Aspect | TK Approximation | IS |
|--------|------------------|-----|
| Speed | Fast | Variable |
| Deterministic | Yes | No |
| Tuning required | Minimal | High |
| High dimensions | Good | Poor |

## References

1. Tierney, L., & Kadane, J. B. (1986). Accurate approximations for posterior moments and marginal densities. *Journal of the American Statistical Association*, 81(393), 82-86.

2. Kass, R. E., & Vos, P. W. (1997). Geometrical foundations of asymptotic inference. *John Wiley & Sons*.

3. Gelman, A., et al. (2013). *Bayesian Data Analysis* (3rd ed.). CRC Press.

## Next Steps

- See "Introduction to TKApprox" for practical usage
- See "Loss Functions" for applying TK to different estimation problems
- See "Simulation Studies" for empirical evaluation of TK accuracy