Je me permets de poser une question car je bloque et je n'arrive pas à avancer sur le sujet.
Je vais partir sur un exemple, monter ce que je souhaite obtenir et où ça bloque. Mon code R utilise le package Weibulltools pour estimer les paramètres de la Weibull et traçage du graphique.
En gros, j'ai 50 avions qui volent depuis deux ans avec une date de début précise (que j'aurais choisi, par exemple, janvier 2016) et une date de fin (décembre 2018), ils ont accumulé 150 000 heures de vol et ont eu un total de 250 pannes sur une pièces lambda (une pompe par exemple).
Ces information me permettent de calculer un MTBF (mean time between failure) qui est égal à
Code : Tout sélectionner
Total des heures de vol / Total des pannes
Dans mon cas, je voudrais obtenir quelques chose de similaire en partant d'une modélisation de Weibull.
J'ai mes données, avec des censures (cas non complet), le MTTF (mean time to failure) est calculable via les paramètres de la Weibull mais pas moyen d'approcher le MTBF. Ma question, est ce possible de le calculer où ca n'a pas de sens ? pas possible ? pour des beta différent de 1 (cas non exponencielle)
Merci bien
CF : http://www.reliawiki.org/index.php/The_Weibull_Distribution
Mon code :
Code : Tout sélectionner
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)