Numerical error in the Weibull model with unexpected posterior value for alpha

Modeling
1

Hi all,

I’m trying to fit a Weibull distribution for some simulated censored time-to-failure data. I get a numerical error error;

Exception thrown at line 26: weibull_log: Shape parameter is 0, but must be > 0! 16

Also, some of my posterior means of the scale are not closer to the fixed values at all (my fixed value of scale is 100 and sometimes I get decimals as the estimated posterior means). Could someone please help me to improve my coding to get the desired output? My attempt of simulation is below.

require (rstan)
require (coda)
require (ggplot2)
require (fitdistrplus)

   model_weibull <- "
    data {
    int<lower=0> N; // number of observations
    int<lower=0> ne; //number of failed observations
    int<lower=0> nc; //number of censored observations
    vector<lower=0> [N]sampledata; //total observations
    vector<lower=0> [ne]failed_obs; //failures
    vector<lower=0> [nc]censored_obs; // Censored observations

}

parameters {
real<lower=0> alpha; // scale
real<lower=0> beta; // shape

}


model {
//define priors of alpha and beta
alpha ~ uniform (0, 150);
beta ~ inv_gamma(0.1, 0.1);

//Likelihood of the data
failed_obs ~ weibull (beta, alpha);
increment_log_prob(weibull_ccdf_log(censored_obs, beta, alpha));
}
"
mod_weibull <- stan_model(model_code = model_weibull)


weibull.data <- function (samplesize, #sample size
                          scale, #scale parameter
                          shape,
                          propFailed){
  data.wei <- sort (rweibull (samplesize, shape = shape, scale = scale))
  T.cen <- data.wei[max(2, ceiling(propFailed*samplesize))]
  censored <- data.wei > T.cen
  event <- c (rep(1, length(data.wei[censored == FALSE])), rep(0, length(data.wei[censored == TRUE])))
  gen.fails <- c (data.wei[censored == FALSE], rep(NA, length(data.wei[censored == TRUE])))
  rank.fails <- rank (gen.fails, ties.method = "first", na.last = "keep")
  gen.fails [is.na(gen.fails)] <- T.cen
  simul.data <- data.frame(sampledata = data.wei,
                           Censored = censored,
                           obs.fails = gen.fails,
                           ranks = rank.fails)
  
  sim.data.bay <- list (sampledata = data.wei,
                        failed_obs = gen.fails[censored==FALSE],
                        censored_obs = gen.fails[censored==TRUE],
                        N = length (data.wei),
                        ne = length(data.wei[censored == FALSE]),
                        nc = length(data.wei[censored == TRUE]))
  
  return(sim.data.bay)
}
num.sim <- 5
results.table <- expand.grid (num.sim = 1:num.sim,
                              samplesize = c (50, 100, 1000),
                              scale.true = c (100),
                              shape.true = c (0.5, 0.75, 1, 1.5, 3),
                              p = c (0.1, 0.8),
                              propFailed = c(0.9))

simulation.bay <-
  function(i) {
    seed = 1000
    set.seed(i * 125 + seed) 
    sim.data <- with(results.table[i,],
                     weibull.data(samplesize = samplesize,
                                  propFailed = propFailed,
                                  shape = shape.true,
                                  scale = scale.true))
    
    
    
    mcmc.weibull <- sampling (mod_weibull, data = sim.data, chains = 1, iter = 10^3,
                              warmup = 50, thin = 5)
    
    
    weibull.coda <- mcmc.list(lapply(1:ncol(mcmc.weibull), 
                                     function(x) {mcmc(as.array(mcmc.weibull)[,x,])}))
    
    posterior.samples <- data.frame (weibull.coda[[1]][, 1:2])
    
    estim.bay <- apply (posterior.samples, 2, median)
    
    Q.true <- with (results.table[i, ],
                    qweibull(p = p,
                             shape = shape.true,
                             scale = scale.true))
    
    
    return (c( estim.bay, Q.true, mcmc.weibull, weibull.coda))
    
  }


sim.results.bay <- lapply(1:nrow(results.table), simulation.bay)
sim.results.bay


scale.bay <- NULL
shape.bay <- NULL
Q.true <- NULL
for (i in 1:nrow(results.table)) {
  
  scale.bay[i] <- (sim.results.bay[[i]]$alpha)
  shape.bay[i] <- (sim.results.bay[[i]][2]$beta)
  Q.true [i] <- sim.results.bay[[i]][[3]]
}


results.table.bay <- cbind (results.table, Q.true, shape.bay, scale.bay)
View (results.table.bay)

Thank you in advance!