The hardware and bandwidth for this mirror is donated by dogado GmbH, the Webhosting and Full Service-Cloud Provider. Check out our Wordpress Tutorial.
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[@]dogado.de.

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.
Health stats visible at Monitor.