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

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

## Introduction

This vignette demonstrates the application of TKApprox to real-world data analysis. We'll analyze a classic reliability dataset using the Weibull distribution with various censoring schemes and loss functions.

## Dataset: Air Conditioning System Failure Times

We'll use the air conditioning system failure times dataset from Proschan (1963), a classic dataset in reliability analysis. The data represents the time intervals between failures of air conditioning systems in aircraft.

```{r}
# Air conditioning failure times (in hours)
ac_failures <- c(23, 261, 87, 7, 120, 14, 62, 47, 225, 71, 246, 21, 42, 20, 5, 
                 12, 120, 11, 3, 14, 71, 11, 14, 11, 16, 90, 1, 16, 52, 95)

cat("Number of observations:", length(ac_failures), "\n")
cat("Mean:", mean(ac_failures), "\n")
cat("Median:", median(ac_failures), "\n")
cat("Range:", range(ac_failures), "\n")
```

## Exploratory Data Analysis

```{r}
# Histogram
hist(ac_failures, breaks = 15, main = "Air Conditioning Failure Times",
     xlab = "Time (hours)", col = "lightblue", freq = FALSE)
lines(density(ac_failures), col = "red", lwd = 2)

# Summary statistics
summary(ac_failures)
```

## Weibull Distribution Model

The Weibull distribution is commonly used for reliability data due to its flexibility in modeling increasing, decreasing, or constant failure rates.

### Define the Distribution

```{r}
# Weibull PDF
pdf_weibull <- function(x, param) {
  dweibull(x, shape = param[1], scale = param[2])
}

# Weibull CDF
cdf_weibull <- function(x, param) {
  pweibull(x, shape = param[1], scale = param[2])
}
```

### Specify Priors

We'll use weakly informative Gamma priors for both parameters:

```{r}
prior_spec <- list(
  shape = list(family = "gamma", hyperparameters = list(shape = 2, rate = 1)),
  scale = list(family = "gamma", hyperparameters = list(shape = 2, rate = 0.01))
)
```

### Fit the Model with Complete Data

```{r}
fit_complete <- tk_fit(
  data = ac_failures,
  censoring_scheme = "complete",
  pdf = pdf_weibull,
  cdf = cdf_weibull,
  prior_spec = prior_spec,
  initial_values = c(shape = 1, scale = 50),
  loss_function = "sel"
)

summary(fit_complete)
```

### Examine Results

```{r}
# Parameter estimates
estimates <- coef(fit_complete)
cat("Shape parameter estimate:", estimates[1], "\n")
cat("Scale parameter estimate:", estimates[2], "\n")

# Covariance matrix
vcov_matrix <- vcov(fit_complete)
cat("\nCovariance matrix:\n")
print(vcov_matrix)

# Standard errors
cat("\nStandard errors:\n")
print(fit_complete$standard_errors)

# Credible intervals
cat("\n95% Credible intervals:\n")
print(fit_complete$credible_intervals)
```

### Model Comparison Statistics

```{r}
print_model_comparison(fit_complete)
```

### Visualization

```{r, fig.width=7, fig.height=6}
# Diagnostic plots
plot(fit_complete, which = 1:4)
```

## Analysis with Different Loss Functions

### LINEX Loss

```{r}
fit_linex <- tk_fit(
  data = ac_failures,
  censoring_scheme = "complete",
  pdf = pdf_weibull,
  cdf = cdf_weibull,
  prior_spec = prior_spec,
  initial_values = c(shape = 1, scale = 50),
  loss_function = "linex",
  loss_params = list(c = 0.1)
)

coef(fit_linex)
```

### General Entropy Loss

```{r}
fit_gel <- tk_fit(
  data = ac_failures,
  censoring_scheme = "complete",
  pdf = pdf_weibull,
  cdf = cdf_weibull,
  prior_spec = prior_spec,
  initial_values = c(shape = 1, scale = 50),
  loss_function = "gel",
  loss_params = list(q = 0.5)
)

coef(fit_gel)
```

### Comparison of Loss Functions

```{r}
comparison <- data.frame(
  Parameter = c("shape", "scale"),
  SEL = coef(fit_complete),
  LINEX = coef(fit_linex),
  GEL = coef(fit_gel)
)

print(comparison)
```

## Analysis with Censored Data

In practice, reliability data often involves censoring. Let's simulate right-censored data from this dataset.

### Create Right-Censored Data

```{r}
# Simulate right censoring at 100 hours
censoring_time <- 100
status <- as.numeric(ac_failures <= censoring_time)

cat("Number of observed failures:", sum(status), "\n")
cat("Number of censored observations:", sum(!status), "\n")
```

### Fit with Right Censoring

```{r, warning=FALSE}
fit_censored <- tk_fit(
  data = ac_failures,
  censoring_scheme = "right-censored",
  pdf = pdf_weibull,
  cdf = cdf_weibull,
  prior_spec = prior_spec,
  initial_values = c(shape = 1, scale = 50),
  loss_function = "sel",
  status = status
)

summary(fit_censored)
```

### Compare Complete vs Censored

```{r}
censoring_comparison <- data.frame(
  Parameter = c("shape", "scale"),
  Complete_Data = coef(fit_complete),
  Right_Censored = coef(fit_censored)
)

print(censoring_comparison)
```

## Reliability Function Estimation

The reliability (survival) function for the Weibull distribution is:

$$R(t) = \exp\left[-\left(\frac{t}{\lambda}\right)^k\right]$$

where $k$ is the shape parameter and $\lambda$ is the scale parameter.

```{r}
# Estimate reliability function
reliability_function <- function(t, shape, scale) {
  exp(-(t / scale)^shape)
}

# Compute reliability at various time points
time_points <- c(10, 20, 50, 100, 200, 500)
est_shape <- coef(fit_complete)[1]
est_scale <- coef(fit_complete)[2]

reliability_estimates <- sapply(time_points, function(t) {
  reliability_function(t, est_shape, est_scale)
})

reliability_table <- data.frame(
  Time = time_points,
  Reliability = reliability_estimates
)

print(reliability_table)

# Plot reliability function
t_seq <- seq(0, 500, length.out = 100)
r_seq <- sapply(t_seq, function(t) reliability_function(t, est_shape, est_scale))

plot(t_seq, r_seq, type = "l", lwd = 2, col = "blue",
     xlab = "Time (hours)", ylab = "Reliability",
     main = "Estimated Reliability Function")
abline(h = 0.5, col = "red", lty = 2)
legend("topright", legend = "50% reliability", col = "red", lty = 2)
```

## Hazard Function Estimation

The hazard function for the Weibull distribution is:

$$h(t) = \frac{k}{\lambda}\left(\frac{t}{\lambda}\right)^{k-1}$$

```{r}
# Estimate hazard function
hazard_function <- function(t, shape, scale) {
  (shape / scale) * (t / scale)^(shape - 1)
}

# Compute hazard at various time points
hazard_estimates <- sapply(time_points, function(t) {
  hazard_function(t, est_shape, est_scale)
})

hazard_table <- data.frame(
  Time = time_points,
  Hazard = hazard_estimates
)

print(hazard_table)

# Plot hazard function
h_seq <- sapply(t_seq, function(t) hazard_function(t, est_shape, est_scale))

plot(t_seq, h_seq, type = "l", lwd = 2, col = "darkgreen",
     xlab = "Time (hours)", ylab = "Hazard Rate",
     main = "Estimated Hazard Function")

# Interpret shape parameter
if (est_shape > 1) {
  cat("\nShape parameter > 1: Increasing failure rate (wear-out)\n")
} else if (est_shape < 1) {
  cat("\nShape parameter < 1: Decreasing failure rate (infant mortality)\n")
} else {
  cat("\nShape parameter = 1: Constant failure rate (exponential)\n")
}
```

## Prior Sensitivity Analysis

Let's examine how sensitive our estimates are to the prior specification.

```{r, fig.width=7, fig.height=6}
sensitivity_shape <- tk_sensitivity(
  fit = fit_complete,
  parameter_name = "shape",
  hyperparameter_name = "shape",
  hyperparameter_values = c(0.5, 1, 2, 5, 10)
)

print(sensitivity_shape)
plot(sensitivity_shape)
```

## Prediction

We can use the fitted model to predict future failure times.

```{r}
# Predict density for new time points
new_times <- c(25, 50, 75, 100, 150)
predicted_density <- predict(fit_complete, newdata = new_times, type = "density")

prediction_table <- data.frame(
  Time = new_times,
  Predicted_Density = predicted_density
)

print(prediction_table)

# Predict survival probability
predicted_survival <- predict(fit_complete, newdata = new_times, type = "survival")

survival_table <- data.frame(
  Time = new_times,
  Survival_Probability = predicted_survival
)

print(survival_table)
```

## Model Diagnostics

```{r}
# Residuals plot
plot(fit_complete, which = 6)

# Check convergence
cat("Convergence code:", fit_complete$convergence, "\n")
cat("Iterations:", fit_complete$iterations, "\n")
cat("Gradient norm:", fit_complete$gradient_norm, "\n")
```

## Comparison with Maximum Likelihood Estimation

Let's compare our Bayesian estimates with frequentist MLE estimates.

```{r}
# MLE using R's built-in function
mle_fit <- MASS::fitdistr(ac_failures, densfun = "weibull")

cat("\n=== MLE Estimates ===\n")
print(mle_fit$estimate)

cat("\n=== Bayesian Estimates (SEL) ===\n")
print(coef(fit_complete))

comparison_mle <- data.frame(
  Parameter = c("shape", "scale"),
  MLE = mle_fit$estimate,
  Bayesian_SEL = coef(fit_complete)
)

print(comparison_mle)
```

## Summary and Interpretation

### Key Findings

1. **Shape parameter**: The estimated shape parameter `r round(coef(fit_complete)[1], 3)` indicates the nature of the failure rate:
   - If > 1: Increasing failure rate (wear-out)
   - If < 1: Decreasing failure rate (infant mortality)
   - If = 1: Constant failure rate (exponential)

2. **Scale parameter**: The estimated scale parameter `r round(coef(fit_complete)[2], 3)` represents the characteristic lifetime.

3. **Reliability**: The reliability function shows the probability of survival beyond time t.

4. **Hazard rate**: The hazard function shows the instantaneous failure rate.

### Practical Implications

- **Maintenance scheduling**: Use the reliability estimates to plan preventive maintenance
- **Warranty analysis**: Estimate probability of failure within warranty period
- **Spares provisioning**: Use hazard rates to determine spare parts inventory
- **System design**: Compare different designs using estimated reliability metrics

## References

Proschan, F. (1963). Theoretical explanation of observed decreasing failure rate. *Technometrics*, 5(3), 375-383.

## Next Steps

- See "Introduction to TKApprox" for basic usage
- See "Censoring Schemes" for handling different censoring patterns
- See "Loss Functions" for alternative estimation methods