Alternative proposed by Varin Sacha using the package weibulltools with delta method (a little bit modified to present only the median and the CI): library(weibulltools) x <- rweibull(100000, 2, 2) data <- reliability_data(x = x, status = 1) ml <- ml_estimation(x = data, distribution = "weibull") conf <- confint_fisher(x = ml, b_lives = 0.5, direction = "x") The results are: > as.numeric(conf[conf$prob == 0.5, c(4, 1, 5)]) [1] 1.657227 1.663275 1.669346 The true median was : > 2*log(2)^(1/2) [1] 1.665109 Let try if the average estimate by Fisher method for 1,000 samples is biased or not: library(weibulltools) library(pbmcapply) k <- pbmclapply(1:1000, FUN=function(i) { x <- rweibull(100000, 2, 2) data <- reliability_data(x = x, status = 1) ml <- ml_estimation(x = data, distribution = "weibull") conf <- confint_fisher(x = ml, b_lives = 0.5, direction = "x") return(as.numeric(conf[conf$prob == 0.5, 1])) }, mc.c...