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

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

library(foocafeReliability)
library(rstan)

Consider the LCD projector failures data considered in the continuous data example:

```{R, include = FALSE, eval = FALSE}

?lcd_projector_failures

We generate data from a Weibull distribution
```{R}
data <- rweibull(n = 25, shape = 2, scale = 4)

Then fit a Weibull model to the data

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

  parameters {
    real<lower=0> shape;
    real<lower=0> scale;
  }

  model {
    shape ~ gamma(2, 2);
    scale ~ inv_gamma(0.1,0.1);

    target += weibull_lpdf(y | shape, scale);
  }
"

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


data_list <- list(N = NROW(data),
                  y = data)

weibull_output <- rstan::sampling(weibull_model, data = data_list, chains = 4, iter = 4000, control = list(adapt_delta = 0.9))

Then we fit an exponential model

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

  parameters {
    real<lower=0> lambda;
  }

  model {
    lambda ~ gamma(0.3, 2);

    target += exponential_lpdf(y | lambda);
  }
"

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


data_list <- list(N = NROW(data),
                  y = data)

exponential_output <- rstan::sampling(exponential_model, data = data_list, chains = 4, iter = 4000, control = list(adapt_delta = 0.9))

stan_hist(weibull_output)
stan_hist(exponential_output)

As expected, the BIC is lower for the

weib <- rstan::extract(weibull_output)
shape <- mean(weib$shape)
scale <- mean(weib$scale)
weibull_BIC <- -2*log(prod(dweibull(data, shape, scale))) + log(NROW(data))*2


expo <- rstan::extract(exponential_output)
lambda <- mean(expo$lambda)
exponential_BIC <- -2*log(prod(dexp(data, lambda))) + log(NROW(data))*1