The hardware and bandwidth for this mirror is donated by METANET, the Webhosting and Full Service-Cloud Provider.
If you wish to report a bug, or if you are interested in having us mirror your free-software or open-source project, please feel free to contact us at mirror[@]metanet.ch.

Real Data Example

Your Name

2026-07-31

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.

# 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")
## Number of observations: 30
cat("Mean:", mean(ac_failures), "\n")
## Mean: 59.6
cat("Median:", median(ac_failures), "\n")
## Median: 22
cat("Range:", range(ac_failures), "\n")
## Range: 1 261

Exploratory Data Analysis

# 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)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##     1.0    12.5    22.0    59.6    83.0   261.0

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

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

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

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)
## 
## === Tierney-Kadane Bayesian Estimation Summary ===
## 
## Model Information:
## -----------------
## Censoring scheme: complete 
## Sample size: 30 
## Number of parameters: 2 
## Loss function: sel 
## 
## Optimization Results:
## --------------------
## Method: nlminb 
## Convergence code: 0 
## Iterations: 11 
## Gradient norm: 0 
## Execution time: 0.33 seconds
## 
## Parameter Estimates:
## --------------------
##  Parameter Posterior_Mode Bayes_Estimate Std_Error  CI_Lower  CI_Upper
##      shape      0.8596845      0.8637463  0.118656  0.631185  1.096308
##      scale     55.9301504     61.2073217 12.394173 36.915190 85.499454
## 
## Model Fit Statistics:
## ---------------------
## Log-posterior at mode: -5.29 
## Log-likelihood at mode: -151.9426 
## Prior contribution: -6.7564 
## 
## Posterior Covariance Matrix:
## ---------------------------
##          shape      scale
## shape 0.014079   0.438637
## scale 0.438637 153.615517

Examine Results

# Parameter estimates
estimates <- coef(fit_complete)
cat("Shape parameter estimate:", estimates[1], "\n")
## Shape parameter estimate: 0.8637463
cat("Scale parameter estimate:", estimates[2], "\n")
## Scale parameter estimate: 61.20732
# Covariance matrix
vcov_matrix <- vcov(fit_complete)
cat("\nCovariance matrix:\n")
## 
## Covariance matrix:
print(vcov_matrix)
##            shape       scale
## shape 0.01407924   0.4386369
## scale 0.43863695 153.6155166
# Standard errors
cat("\nStandard errors:\n")
## 
## Standard errors:
print(fit_complete$standard_errors)
##     shape     scale 
##  0.118656 12.394173
# Credible intervals
cat("\n95% Credible intervals:\n")
## 
## 95% Credible intervals:
print(fit_complete$credible_intervals)
##           2.5 %    97.5 %
## shape  0.631185  1.096308
## scale 36.915190 85.499454

Model Comparison Statistics

print_model_comparison(fit_complete)
## 
## === Model Comparison Statistics ===
## 
##                Statistic       Value
##           Log-Likelihood -151.942618
##  Negative Log-Likelihood  151.942618
##                      AIC  307.885236
##                      BIC  310.687631
##                     CAIC  312.687631
##                     HQIC  308.781746
##                      DIC  307.885236
##   Expected Log-Posterior   -5.289968
##     Number of Parameters    2.000000
##              Sample Size   30.000000

Visualization

# Diagnostic plots
plot(fit_complete, which = 1:4)

Analysis with Different Loss Functions

LINEX Loss

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)
##      shape      scale 
##  0.8488107 53.4599974

General Entropy Loss

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)
##      shape      scale 
##  0.8513032 58.8918032

Comparison of Loss Functions

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

print(comparison)
##       Parameter        SEL      LINEX        GEL
## shape     shape  0.8637463  0.8488107  0.8513032
## scale     scale 61.2073217 53.4599974 58.8918032

Analysis with Censored Data

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

Create Right-Censored Data

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

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

Fit with Right Censoring

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)
## 
## === Tierney-Kadane Bayesian Estimation Summary ===
## 
## Model Information:
## -----------------
## Censoring scheme: right-censored 
## Sample size: 30 
## Number of parameters: 2 
## Loss function: sel 
## 
## Optimization Results:
## --------------------
## Method: nlminb 
## Convergence code: 0 
## Iterations: 11 
## Gradient norm: 0 
## Execution time: 0.8488 seconds
## 
## Parameter Estimates:
## --------------------
##  Parameter Posterior_Mode Bayes_Estimate Std_Error   CI_Lower    CI_Upper
##      shape      0.7210545      0.7204095   0.11390  0.4971696   0.9436493
##      scale     66.3057290     76.7894829  18.18255 41.1523366 112.4266292
## 
## Model Fit Statistics:
## ---------------------
## Log-posterior at mode: -4.5292 
## Log-likelihood at mode: -129.1484 
## Prior contribution: -6.7272 
## 
## Posterior Covariance Matrix:
## ---------------------------
##          shape      scale
## shape 0.012973   0.195768
## scale 0.195768 330.605184

Compare Complete vs Censored

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

print(censoring_comparison)
##       Parameter Complete_Data Right_Censored
## shape     shape     0.8637463      0.7204095
## scale     scale    61.2073217     76.7894829

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.

# 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)
##   Time Reliability
## 1   10 0.811295575
## 2   20 0.683483419
## 3   50 0.431831629
## 4  100 0.216950507
## 5  200 0.061992226
## 6  500 0.002163706
# 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}\]

# 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)
##   Time     Hazard
## 1   10 0.01806291
## 2   20 0.01643506
## 3   50 0.01450609
## 4  100 0.01319879
## 5  200 0.01200930
## 6  500 0.01059978
# 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")
}
## 
## Shape parameter < 1: Decreasing failure rate (infant mortality)

Prior Sensitivity Analysis

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

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)
## Prior Sensitivity Analysis
## ==========================
## Parameter: shape 
## Hyperparameter: shape 
## Loss function: sel 
## Number of hyperparameter values tested: 5 
## 
## Results:
##  hyperparameter_value log_posterior log_likelihood convergence iterations
##                   0.5             0              0           0          0
##                   1.0             0              0           0          0
##                   2.0             0              0           0          0
##                   5.0             0              0           0          0
##                  10.0             0              0           0          0
##  estimate_shape se_shape estimate_scale se_scale
##              NA       NA             NA       NA
##              NA       NA             NA       NA
##              NA       NA             NA       NA
##              NA       NA             NA       NA
##              NA       NA             NA       NA
plot(sensitivity_shape)
## True parameter not available; cannot compute risk.

Prediction

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

# 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)
##   Time Predicted_Density
## 1   25       0.010049942
## 2   50       0.006264190
## 3   75       0.004168028
## 4  100       0.002863484
## 5  150       0.001427531
# 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)
##   Time Survival_Probability
## 1   25            0.6303716
## 2   50            0.4318316
## 3   75            0.3036500
## 4  100            0.2169505
## 5  150            0.1142996

Model Diagnostics

# Residuals plot
plot(fit_complete, which = 6)

# Check convergence
cat("Convergence code:", fit_complete$convergence, "\n")
## Convergence code: 0
cat("Iterations:", fit_complete$iterations, "\n")
## Iterations: 11
cat("Gradient norm:", fit_complete$gradient_norm, "\n")
## Gradient norm: 1.322962e-10

Comparison with Maximum Likelihood Estimation

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

# MLE using R's built-in function
mle_fit <- MASS::fitdistr(ac_failures, densfun = "weibull")
## Warning in densfun(x, parm[1], parm[2], ...): NaNs produced
cat("\n=== MLE Estimates ===\n")
## 
## === MLE Estimates ===
print(mle_fit$estimate)
##      shape      scale 
##  0.8530424 54.4732388
cat("\n=== Bayesian Estimates (SEL) ===\n")
## 
## === Bayesian Estimates (SEL) ===
print(coef(fit_complete))
##      shape      scale 
##  0.8637463 61.2073217
comparison_mle <- data.frame(
  Parameter = c("shape", "scale"),
  MLE = mle_fit$estimate,
  Bayesian_SEL = coef(fit_complete)
)

print(comparison_mle)
##       Parameter        MLE Bayesian_SEL
## shape     shape  0.8530424    0.8637463
## scale     scale 54.4732388   61.2073217

Summary and Interpretation

Key Findings

  1. Shape parameter: The estimated shape parameter 0.864 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 61.207 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

References

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

Next Steps

These binaries (installable software) and packages are in development.
They may not be fully stable and should be used with caution. We make no claims about them.