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

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

## Introduction

TKApprox supports multiple Bayesian loss functions for parameter estimation. This vignette explains each loss function, its properties, and when to use it.

## Squared Error Loss (SEL)

Squared Error Loss is the most common loss function, leading to the posterior mean as the Bayes estimator.

**Loss function:** $L(\theta, \hat{\theta}) = (\theta - \hat{\theta})^2$

**Bayes estimator:** $\hat{\theta} = E[\theta | x]$ (posterior mean)

**Properties:**
- Symmetric: overestimation and underestimation are penalized equally
- Unbiased under regularity conditions
- Most commonly used loss function

```{r}
# Define exponential distribution
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 with SEL (default)
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)
```

## LINEX Loss

LINEX (Linear-Exponential) Loss is an asymmetric loss function useful when overestimation and underestimation have different consequences.

**Loss function:** $L(\theta, \hat{\theta}) = [\exp(c(\hat{\theta} - \theta)) - c(\hat{\theta} - \theta) - 1]$

**Bayes estimator:** $\hat{\theta} = -\frac{1}{c} \log E[\exp(-c\theta) | x]$

**Properties:**
- Asymmetric: $c > 0$ penalizes overestimation more, $c < 0$ penalizes underestimation more
- Approaches SEL as $c \to 0$
- Useful in reliability and risk assessment

```{r}
# Fit with LINEX loss (c = 0.5, penalizes overestimation)
fit_linex_pos <- 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)
)

# Fit with LINEX loss (c = -0.5, penalizes underestimation)
fit_linex_neg <- 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)
)

# Compare estimates
data.frame(
  SEL = coef(fit_sel),
  LINEX_c_0.5 = coef(fit_linex_pos),
  LINEX_c_neg0.5 = coef(fit_linex_neg)
)
```

### LINEX Sensitivity Analysis

```{r}
# Examine sensitivity to LINEX parameter c
c_values <- c(-2, -1, -0.5, -0.1, 0.1, 0.5, 1, 2)
linex_estimates <- sapply(c_values, function(c) {
  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 = "linex",
    loss_params = list(c = c)
  )
  coef(fit)
})

plot(c_values, linex_estimates, type = "b", pch = 19,
     xlab = "LINEX parameter c", ylab = "Estimate",
     main = "LINEX Estimates vs c")
abline(h = coef(fit_sel), col = "red", lty = 2)
legend("topright", legend = c("LINEX", "SEL"), col = c("black", "red"),
       pch = c(19, NA), lty = c(1, 2))
```

## General Entropy Loss (GEL)

General Entropy Loss is another asymmetric loss function useful for scale parameters.

**Loss function:** $L(\theta, \hat{\theta}) = (\hat{\theta}/\theta)^q / q - \log(\hat{\theta}/\theta) - 1/q$

**Bayes estimator:** $\hat{\theta} = [E(\theta^{-q} | x)]^{-1/q}$

**Properties:**
- Asymmetric: $q > 0$ penalizes overestimation, $q < 0$ penalizes underestimation
- Scale-invariant
- Useful for positive parameters (rates, scales)

```{r}
# Fit with GEL (q = 0.5)
fit_gel_pos <- 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)
)

# Fit with GEL (q = -0.5)
fit_gel_neg <- 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)
)

# Compare estimates
data.frame(
  SEL = coef(fit_sel),
  GEL_q_0.5 = coef(fit_gel_pos),
  GEL_q_neg0.5 = coef(fit_gel_neg)
)
```

## Precautionary Loss

Precautionary Loss is useful when we want to be conservative in estimation.

**Loss function:** $L(\theta, \hat{\theta}) = \hat{\theta}^2/\theta - 2\hat{\theta} + \theta$

**Bayes estimator:** $\hat{\theta} = \sqrt{E(\theta^2 | x)}$

**Properties:**
- Conservative: tends to produce larger estimates
- Useful for safety-critical applications
- Only appropriate for positive parameters

```{r}
# Fit with precautionary loss
fit_precautionary <- tk_fit(
  data = data,
  censoring_scheme = "complete",
  pdf = pdf_exp,
  cdf = cdf_exp,
  prior_spec = prior_spec,
  initial_values = c(rate = 1),
  loss_function = "precautionary"
)

# Compare with SEL
data.frame(
  SEL = coef(fit_sel),
  Precautionary = coef(fit_precautionary)
)
```

## Weighted Squared Error Loss

Weighted SEL uses inverse-variance weighting, useful when parameters have different scales.

**Loss function:** $L(\theta, \hat{\theta}) = (\theta - \hat{\theta})^2 / \theta^2$

**Bayes estimator:** $\hat{\theta} = 1 / E(\theta^{-1} | x)$

**Properties:**
- Scale-invariant
- Useful for parameters with large variability
- Emphasizes relative error

```{r}
# Fit with weighted SEL
fit_weighted <- tk_fit(
  data = data,
  censoring_scheme = "complete",
  pdf = pdf_exp,
  cdf = cdf_exp,
  prior_spec = prior_spec,
  initial_values = c(rate = 1),
  loss_function = "weighted-sel"
)

# Compare with SEL
data.frame(
  SEL = coef(fit_sel),
  Weighted_SEL = coef(fit_weighted)
)
```

## Custom Loss Functions

You can also specify custom loss functions by providing the g(θ) function directly.

```{r}
# Example: Estimate the median of the posterior
# For exponential distribution, median = log(2)/rate
# We want to estimate log(rate) instead of rate directly
custom_g <- function(param) log(param[1])

fit_custom <- tk_fit(
  data = data,
  censoring_scheme = "complete",
  pdf = pdf_exp,
  cdf = cdf_exp,
  prior_spec = prior_spec,
  initial_values = c(rate = 1),
  loss_function = "custom",
  custom_g = custom_g
)

# The estimate is E[log(rate) | x]
# Transform back to rate scale
rate_estimate <- exp(coef(fit_custom))

data.frame(
  SEL_rate = coef(fit_sel),
  Custom_log_rate = rate_estimate
)
```

## Comparing All Loss Functions

```{r}
# Fit with all loss functions
fit_sel <- tk_fit(data, "complete", pdf_exp, cdf_exp, prior_spec,
                  initial_values = c(rate = 1), loss_function = "sel")

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

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

fit_precautionary <- tk_fit(data, "complete", pdf_exp, cdf_exp, prior_spec,
                            initial_values = c(rate = 1), loss_function = "precautionary")

fit_weighted <- tk_fit(data, "complete", pdf_exp, cdf_exp, prior_spec,
                       initial_values = c(rate = 1), loss_function = "weighted-sel")

# Compare all estimates
comparison <- data.frame(
  Loss_Function = c("SEL", "LINEX (c=0.5)", "GEL (q=0.5)", "Precautionary", "Weighted SEL"),
  Estimate = c(coef(fit_sel), coef(fit_linex), coef(fit_gel),
               coef(fit_precautionary), coef(fit_weighted))
)

print(comparison)

# Visual comparison
barplot(comparison$Estimate, names.arg = comparison$Loss_Function,
        main = "Bayes Estimates Under Different Loss Functions",
        ylab = "Rate Estimate", col = "steelblue")
abline(h = 1.5, col = "red", lty = 2)  # True value
legend("topright", legend = "True value", col = "red", lty = 2)
```

## Choosing a Loss Function

The choice of loss function depends on your application:

- **SEL**: Use when errors are symmetric and you want unbiased estimates
- **LINEX**: Use when overestimation and underestimation have different costs
  - $c > 0$: Overestimation is more costly (e.g., safety margins)
  - $c < 0$: Underestimation is more costly (e.g., resource allocation)
- **GEL**: Use for scale parameters where relative error matters
- **Precautionary**: Use when you want conservative estimates
- **Weighted SEL**: Use when parameters have different scales or variability

## Multi-Parameter Models

For multi-parameter models, each parameter can be estimated under a different loss function:

```{r}
# Two-parameter Weibull example
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 = "gamma", hyperparameters = list(shape = 2, rate = 1)),
  scale = list(family = "gamma", hyperparameters = list(shape = 2, rate = 1))
)

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

# Fit with SEL (applies to both parameters)
fit_weibull <- tk_fit(
  data = data_weibull,
  censoring_scheme = "complete",
  pdf = pdf_weibull,
  cdf = cdf_weibull,
  prior_spec = prior_spec,
  initial_values = c(shape = 1.5, scale = 1),
  loss_function = "sel"
)

coef(fit_weibull)
```

## Next Steps

- See "Prior Specification" for information on choosing priors
- See "Simulation Studies" for comparing loss functions in simulation