---
title: "Prior Specification"
author: "Your Name"
date: "`r Sys.Date()`"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Prior Specification}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

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

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

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

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

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

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

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

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

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

### Exponential Prior

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

**Parameters:** rate (λ)

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

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

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

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

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

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

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

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

## Independent Priors for Multiple Parameters

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

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

## Custom Prior Functions

You can also specify a custom prior function directly:

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

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

## Prior Sensitivity Analysis

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

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

## Systematic Prior Sensitivity

Use `tk_sensitivity()` for systematic examination of prior hyperparameters:

```{r}
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)
plot(sensitivity)
```

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

1. **Check prior support**: Ensure the prior support matches the parameter domain
2. **Scale matters**: Consider the scale of your parameters when setting hyperparameters
3. **Sensitivity analysis**: Always check how sensitive results are to prior choices
4. **Conjugate when possible**: Use conjugate priors for computational efficiency
5. **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