Теория Вейбулла с расчетом R: MTBF для данных цензоров

Question

Мне нужна помощь по вопросу в теории Вейбулла, я приведу пример, чтобы проиллюстрировать это. Предположим, что у нас есть 50 самолетов, которые летят с 2016 по декабрь 2018 года (допустим, общее количество летных часов составляет 150 000), и на этом самолете было 250 съемок конкретного насоса. Основным методом, если я хочу рассчитать MTBF (среднее время между отказами) насоса для всего парка, я делаю:total flight hours / total removalsв моем случае это будет 150 000/250 = 600 hours, Теперь, что мне нужно, это тот же результат из оценки Вейбулла. Мне нужно что-то вроде формулы MTTF (среднее время до отказа) для MTBF в случае, когда у меня есть данные цензоров. Например,

list_of_packages <- c("weibulltools", "lubridate", "dplyr")

new_packages <- list_of_packages[!(list_of_packages %in% installed.packages()[,"Package"])]

if(length(new_packages)){
  install.packages(new_packages)
} 

lapply(list_of_packages, require, character.only = TRUE)

#### Simulation weibull vs cumulative MTBF ####

## cass Beta > 1
set.seed(1234)
n <- 1000
b <- 5
e <- 5000

x <- rweibull(n, b, e)

## No censor data, mttf and MTBF
mttf_weib <- round(e * gamma(1 + (1 / b)), 0)
mtbf_cum <- sum(x)/n

mttf_weib
mtbf_cum 


## With Censor data
y <- rweibull(n, b, e)

df <- data.frame(x = x, y = y) %>% 
  mutate(time = pmin(x, y),
         event = as.numeric(x < y),
         ac = as.character("ouistiti"))

data <- df
# Weibulltools package 
distribution <- "weibull"

df_meth <- johnson_method(x = data$time, event = data$event, id = data$ac)
df_meth


df_est <- ml_estimation(x = df_meth$characteristic,
                        event = df_meth$status,
                        distribution = distribution, 
                        conf_level = 0.90)

df_est <- rank_regression(x = df_meth$characteristic,
                          y = df_meth$prob,
                          event = df_meth$status,
                          distribution = distribution,
                          conf_level = .95)


conf_betabin <- confint_betabinom(x = df_meth$characteristic,
                                  event = df_meth$status,
                                  loc_sc_params = df_est$loc_sc_coefficients,
                                  distribution = distribution,
                                  bounds = "two_sided",
                                  conf_level = 0.95,
                                  direction = "y")

plot_weibull <- plot_prob(x = df_meth$characteristic,
                          y = df_meth$prob,
                          event = df_meth$status,
                          id = df_meth$id,
                          distribution = distribution,
                          title_main = "Weibull Analysis",
                          title_x = "Flight hours",
                          title_y = "Probability of Failure",
                          title_trace = "Removed Items")

plot_reg_weibull <- plot_mod(p_obj = plot_weibull,
                             x = conf_betabin$characteristic,
                             y = conf_betabin$prob,
                             loc_sc_params = df_est$loc_sc_coefficients,
                             distribution = distribution,
                             title_trace = "Estimated Weibull CDF")

plot_conf_beta <- plot_conf(p_obj = plot_reg_weibull,
                            x = list(conf_betabin$characteristic),
                            y = list(conf_betabin$lower_bound,
                                     conf_betabin$upper_bound),
                            direction = "y",
                            distribution = distribution,
                            title_trace = "Confidence Region")
#### Suppression warnings
plot_conf_beta$elementId <- NULL

ls <- list(plot_weibull = plot_conf_beta, 
           n_failure = paste("Nombre de failure" , nrow(data)),
           beta = df_est$coefficients[2],
           eta = df_est$coefficients[1],
           MTTF = round(round(as.numeric(df_est$coefficients[1]), 2) * gamma(1 + (1 / round(as.numeric(df_est$coefficients[2]), 2))), 0),
           data = data)
ls

#Proportion of censors data
table(df$event)
# beta
ls$beta
# eta
ls$eta
# MTBF
ls$MTTF

# 2 times greather ??
sum(df$time)/sum(df$event)

Я понимаю, что, возможно, это не очень хорошая тема, но я не получил ответов на перекрестном форуме. Спасибо