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

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

## Introduction

TKApprox provides a distribution-independent framework for Bayesian estimation of arbitrary univariate probability models using the Tierney-Kadane approximation. This vignette provides a gentle introduction to the package's main features and workflow.

## Basic Workflow

The typical workflow in TKApprox follows these steps:

1. Define your probability distribution (PDF/PMF and CDF)
2. Specify prior distributions for parameters
3. Provide your data and censoring scheme
4. Call `tk_fit()` to perform Bayesian estimation
5. Examine results using S3 methods and visualization

## A Simple Example: Exponential Distribution

Let's start with a simple example using the exponential distribution.

### Step 1: Define the Distribution

```{r}
# Define the exponential PDF
pdf_exp <- function(x, param) {
  dexp(x, rate = param)
}

# Define the exponential CDF
cdf_exp <- function(x, param) {
  pexp(x, rate = param)
}
```

### Step 2: Specify Prior

```{r}
# Gamma prior for the rate parameter
prior_spec <- list(
  rate = list(family = "gamma", hyperparameters = list(shape = 2, rate = 1))
)
```

### Step 3: Generate Data

```{r}
set.seed(123)
data <- rexp(20, rate = 1.5)
```

### Step 4: Fit the Model

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

### Step 5: Examine Results

```{r}
# Print basic results
print(fit)

# Detailed summary
summary(fit)

# Extract Bayes estimates
coef(fit)

# Extract covariance matrix
vcov(fit)

# Model comparison statistics
print_model_comparison(fit)
```

## Visualization

```{r}
# Plot all diagnostics
plot(fit)

# Or select specific plots
plot(fit, which = 1)  # Posterior approximation
plot(fit, which = 2)  # Likelihood surface
plot(fit, which = 3)  # Prior vs posterior
```

## Different Loss Functions

TKApprox supports multiple Bayesian loss functions:

### Squared Error Loss (Posterior Mean)

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

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

### General Entropy Loss

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

## Censored Data

TKApprox handles various censoring schemes. Here's an example with right-censored data:

```{r}
# Create right-censored data
status <- c(1, 1, 0, 1, 0, 1, 1, 0, 1, 1)  # 1 = observed, 0 = censored

fit_censored <- tk_fit(
  data = data,
  censoring_scheme = "right-censored",
  pdf = pdf_exp,
  cdf = cdf_exp,
  prior_spec = prior_spec,
  initial_values = c(rate = 1),
  loss_function = "sel",
  status = status
)

summary(fit_censored)
```

## Prior Sensitivity Analysis

Examine how sensitive your estimates are to prior hyperparameters:

```{r}
sensitivity <- tk_sensitivity(
  fit = fit,
  parameter_name = "rate",
  hyperparameter_name = "shape",
  hyperparameter_values = c(0.5, 1, 2, 5, 10)
)

print(sensitivity)
plot(sensitivity)
```

## Multi-Parameter Models

TKApprox works with models of any dimensionality. Here's a two-parameter example with the Weibull distribution:

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

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

# Fit the model
fit_weibull <- tk_fit(
  data = data_weibull,
  censoring_scheme = "complete",
  pdf = pdf_weibull,
  cdf = cdf_weibull,
  prior_spec = prior_spec_weibull,
  initial_values = c(shape = 1.5, scale = 1),
  loss_function = "sel"
)

summary(fit_weibull)
plot(fit_weibull, which = 2)  # Likelihood surface for 2-parameter model
```

## Next Steps

- See "Defining Custom Distributions" for more complex distribution examples
- See "Censoring Schemes" for detailed coverage of all censoring types
- See "Loss Functions" for more information on Bayesian loss functions
- See "Prior Specification" for advanced prior modeling