In a-segar/foocafeReliability: Bayesian Reliability Tools and Examples

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

library(foocafeReliability)
library(rstan)

Ku et al. (1972) reports on fatigue testing of bearings used with a particular lubricant and assumes that the failure times follow a Weibull distribution. The experimenters used 10 testers, bench-type rigs, and found that the testers impacted the measured failure times. The dataset 'lubricant_bearing_failures' presents the bearing failure time data (in hours) that they collected when they used an aviation gas turbine lubricant O-64-2. The experimenters want to determine the bearing failure time distribution when they use the bearings with O-64-2, by removing the tester effect. In analysing these data, we can account for the tester-to-tester differences by specifying

$$ Y \sim Weibull(\alpha_i, \beta) $$

where $y_{ij}$ is the $j-th$ observed failure time from the $i-th$ tester. Note that this model specification used the first parametrisation of the Weibull distribution given in Appendix B. Further, we model the logged scale parameter $\alpha_i$ by:

lubricant_bearing_failures

Used example from https://www.r-bloggers.com/hierarchical-models-with-rstan-part-1/

Model from book (couldn't get the scale parameter results to match the books based on Ex. 4.8 Bayesian Relaibility):

The goal is to show how to model with different subgroups

modelString = "
  data {
    int N; //the number of observations
    int J; //the number of testers
    int id[N]; //vector of tester indices
    vector[N] y; //the response variable
  }

  parameters {
    real<lower=0> sigma_2;

    vector[J] gamma;
    real mu;

    real<lower=0> beta;
  }

  transformed parameters {
    real sigma;

    sigma = sqrt(sigma_2);
  }

  model {
    vector[N] log_alpha;
    vector[N] alpha;
    vector[N] scale;

    //priors
    sigma_2 ~ inv_gamma(0.01, 0.01);

    mu ~ normal(0, sqrt(1000));
    beta ~ gamma(1.5, 0.5);

    for(j in 1:J){
      gamma[j] ~ normal(0, sigma); 
    }

    for(n in 1:N){
      log_alpha[n] = mu + gamma[id[n]];
      alpha[n] = exp(log_alpha[n]);

      scale[n] = pow(alpha[n], beta);
    }

    //likelihood
    target += weibull_lpdf(y | beta, scale);
  }
"

lubricant_bearing_failures_model <- rstan::stan_model(model_code = modelString)


lubricant_bearing_fail_list <- list(N = NROW(lubricant_bearing_failures$failure_times),
                            J = NROW(unique(lubricant_bearing_failures$tester)),
                            id = lubricant_bearing_failures$tester,
                            y = lubricant_bearing_failures$failure_times
                          )

output <- rstan::sampling(lubricant_bearing_failures_model, data = lubricant_bearing_fail_list, chains = 4, iter = 4000, control = list(adapt_delta = 0.9))

Experimenting with priors

modelString = "
  data {
    int N; //the number of observations
    int J; //the number of testers
    int id[N]; //vector of tester indices
    vector[N] y; //the response variable
  }

  parameters {
    real<lower=0> test_prior;
    real<lower=0> test_prior_lnorm;
    real<lower=0> test_prior_gamma;
    real<lower=0> test_prior_gamma2;
    real<lower=0> test_prior_uniform;
  }

  model {
    test_prior ~ gamma(2.5,2350);
    test_prior_lnorm ~ lognormal(5, 0.5);
    test_prior_gamma ~ gamma(2.5, 0.6);
    test_prior_gamma2 ~ gamma(2, 1);
    test_prior_uniform ~ uniform(0,100);
  }
"
lubricant_bearing_failures_model2 <- rstan::stan_model(model_code = modelString)


lubricant_bearing_fail_list <- list(N = NROW(lubricant_bearing_failures$failure_times),
                            J = NROW(unique(lubricant_bearing_failures$tester)),
                            id = lubricant_bearing_failures$tester,
                            y = lubricant_bearing_failures$failure_times
                          )

output2 <- rstan::sampling(lubricant_bearing_failures_model2, data = lubricant_bearing_fail_list, chains = 4, iter = 4000, control = list(adapt_delta = 0.9))

Model 1

library(foocafeReliability)

modelString = "
  data {
    int N; //the number of observations
    int J; //the number of testers
    int id[N]; //vector of tester indices
    vector[N] y; //the response variable
  }

  parameters {
    vector<lower=0>[J] alpha;
    real<lower=0> beta;

    real mu;
    real sigma;
  }

  model {
    vector[N] scale;

    beta ~ gamma(1.5, 0.5);
    mu ~ gamma(2,1);
    sigma ~ uniform(0, 100);

    for(j in 1:J){
      alpha[j] ~ lognormal(mu, sigma);
    }

    for (n in 1:N)
      scale[n] = alpha[id[n]];

    //likelihood
    target += weibull_lpdf(y | beta, scale);
  }
"

lubricant_bearing_failures_model <- rstan::stan_model(model_code = modelString)


lubricant_bearing_fail_list <- list(N = NROW(lubricant_bearing_failures$failure_times),
                            J = NROW(unique(lubricant_bearing_failures$tester)),
                            id = lubricant_bearing_failures$tester,
                            y = lubricant_bearing_failures$failure_times
                          )

output <- rstan::sampling(lubricant_bearing_failures_model, data = lubricant_bearing_fail_list, chains = 4, iter = 4000, control = list(adapt_delta = 0.9))

supercomputer_failures_model

Model 2

modelString = "
  data {
    int N; //the number of observations
    int J; //the number of testers
    int id[N]; //vector of tester indices
    vector[N] y; //the response variable
  }

  parameters {
    vector<lower=0>[J] alpha;
    real<lower=0> beta;

    real mu;
    real sigma;
  }

  transformed parameters {
    vector[N] scale;

    for (n in 1:N)
      scale[n] = alpha[id[n]];
  }

  model {
    beta ~ gamma(1.5, 0.5);
    mu ~ gamma(2,1);
    sigma ~ inv_gamma(0.1, 0.1);

    for(j in 1:J){
      alpha[j] ~ lognormal(mu, sigma);
    }

    //likelihood
    target += weibull_lpdf(y | beta, scale);
  }

  generated quantities {
    //sample predicted values from the model for posterior predictive checks
    real y_rep[N];

    for(n in 1:N)
      y_rep[n] = weibull_rng(beta, scale[n]);
  }
"

lubricant_bearing_failures_model2 <- rstan::stan_model(model_code = modelString)


lubricant_bearing_fail_list <- list(N = NROW(lubricant_bearing_failures$failure_times),
                            J = NROW(unique(lubricant_bearing_failures$tester)),
                            id = lubricant_bearing_failures$tester,
                            y = lubricant_bearing_failures$failure_times
                          )

output2 <- rstan::sampling(lubricant_bearing_failures_model2, data = lubricant_bearing_fail_list, chains = 4, iter = 4000, control = list(adapt_delta = 0.9))

costTableList <- list()
dataList <- list()

t = c(50, 100, 150, 200)

output_vals <- rstan::extract(output)

for (kk in 1:10){
  for (jj in 1:NROW(output_vals)){

  thisShape <- output_vals$beta[jj]
  thisScale <- output_vals$alpha[kk,jj]

  get_reliability <- function(thisShape, thisScale){
    f1 <- function(t){
      exp(-((t/thisScale)^thisShape))
    }
    return(f1)
  }

  reliability <- get_reliability(thisShape, thisScale)

  costTable <- data.frame("t" = t,
                          "reliability" = numeric(NROW(t)))

  costTable$reliability <- reliability(t)

  costTableList[[jj]] <- costTable
  }

  temp <- lapply(costTableList, '[', 2) %>% dplyr::bind_cols()

  data <- data.frame(t = t,
                     meanVal = apply(temp, 1, mean),
                     LQ = apply(temp, 1, function(x) {quantile(x, 0.05)}),
                     UQ = apply(temp, 1, function(x) {quantile(x, 0.95)}))

  dataList[[kk]] <- data

}

data <- plyr::aaply(plyr::laply(dataList, as.matrix), c(2, 3), mean) %>%
  as.data.frame()

ggplot(data = data, aes(t)) +
  geom_line(aes(y = meanVal)) +
  geom_ribbon(aes(ymin = LQ, ymax = UQ), alpha = 0.05)

Evaluate goodness of fit

n <- NROW(lubricant_bearing_failures)
K <- round(n^0.4)
a <- seq(0,1,1/K)
p <- diff(a)

probabilities <- numeric(n)
m <- numeric(K)
m_temp <- numeric()
output_vals <- rstan::extract(output)

chi_val <- qchisq(0.95,K-1)
r_b <- numeric(NROW(output_vals[[1]]))

for (jj in 1:NROW(r_b)) {

  thisShape <- output_vals$beta[jj]
  thisScale <- output_vals$alpha[jj,]

  for (ii in 1:length(lubricant_bearing_failures$failure_times)){
    probabilities[ii] <- pweibull(lubricant_bearing_failures$failure_times[ii],
                                  thisShape,
                                  thisScale[lubricant_bearing_failures$tester[ii]])
  }

  for (ii in 1:length(m)){
    m[ii] <- sum(probabilities > a[ii] & probabilities < a[ii+1])
  }

  r_b[jj] <- sum(((m - n*p)^2)/(n*p))

  m_temp <- rbind(m_temp, m)

}

sum(r_b > chi_val)/NROW(r_b) * 100