## ----setup, include=FALSE-----------------------------------------------------
knitr::opts_chunk$set(echo = TRUE)
library(TKApprox)

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

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

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

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

## -----------------------------------------------------------------------------
print_model_comparison(fit_complete)

## ----fig.width=7, fig.height=6------------------------------------------------
# Diagnostic plots
plot(fit_complete, which = 1:4)

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

## -----------------------------------------------------------------------------
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 <- data.frame(
  Parameter = c("shape", "scale"),
  SEL = coef(fit_complete),
  LINEX = coef(fit_linex),
  GEL = coef(fit_gel)
)

print(comparison)

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

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

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

print(censoring_comparison)

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

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

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

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

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

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