请求协助：将Birnbaum-Saunders分布拟合至右删失数据

需要用R语言对客户等待时间的右删失数据拟合Birnbaum-Saunders分布，数据规则为：0表示客户未得到服务，等待时间未删失；1表示客户已得到服务，等待时间为右删失。已尝试fitdistcens、flexsurvreg等多个R包与函数，但未能实现目标。以下是部分尝试代码：

尝试1 #

g1 <- read_excel("/Users/angelavidalmonge/Desktop/TFG/Datos/datos_censurados_7_9h.xlsx", 
                 sheet = 1, col_names = FALSE)
colnames(g1)[1] <- "TiemposEspera"
colnames(g1)[2] <- "Censura"
cens <- g1$Censura == 1

library(dplyr)

g1_cens_interv <- g1 %>%
  mutate(
    left = ifelse(Censura == 0, TiemposEspera, TiemposEspera),   # Ambos igual
    right = ifelse(Censura == 0, TiemposEspera, NA)             # Si censurado, intervalo abierto
  ) %>%
  select(left, right)
g1_cens_interv <- as.data.frame(g1_cens_interv)

# Crear una función lista para fitdistcens con nombre "bisa"
dbs  <- function(x, shape, scale) dbisa(x, shape = shape, scale = scale)
pbs   <- function(q, shape, scale) pbisa(q, shape = shape, scale = scale)
dbs(1, shape = 1, scale = 1)   # debería devolver un valor numérico
pbs(1, shape = 1, scale = 1)

ajuste_bs <- fitdistcens(
  censdata = g1_cens_interv,
  distr = "bs",
  start = list(shape = 1, scale = mean(g1_cens_interv$left))
)
shape2 <- ajuste_bs$estimate["shape"]
scale2 <- ajuste_bs$estimate["scale"]

尝试2 #

library(VGAM)

dBS <- function(x, alpha, beta, log = FALSE) {
  dbisa(x, shape = alpha, scale = beta, log = log)
}

# Función de distribución acumulativa (CDF)
pBS <- function(q, alpha, beta, lower.tail = TRUE, log.p = FALSE) {
  pbisa(q, shape = alpha, scale = beta, lower.tail = lower.tail, log.p = log.p)
}

# Función de cuantiles (inversa de CDF)
qBS <- function(p, alpha, beta, lower.tail = TRUE, log.p = FALSE) {
  qbisa(p, shape = alpha, scale = beta, lower.tail = lower.tail, log.p = log.p)
}

# Generación de números aleatorios
rBS <- function(n, alpha, beta) {
  rbisa(n, shape = alpha, scale = beta)
}

birnbaum_saunders <- list(
  name = "BS",
  pars = c("alpha", "beta"),
  location = NULL,  # parámetro principal (escala)
  transforms = c(log, log),
  inv.transforms = c(exp, exp),
  inits = function(t) {
    c(1, mean(t))  # valores iniciales razonables
  }
)
?flexsurvreg
ajuste_bs <- flexsurvreg(g1_cens ~ 1, data = g1, dist = birnbaum_saunders)
alpha_bs <- as.numeric(exp(ajuste_bs$coefficients["alpha"]))
beta_bs <- as.numeric(exp(ajuste_bs$coefficients["beta"]))
alpha_bs
beta_bs

若需要，可提供数据Excel文件。此内容为重要毕业设计内容，恳请协助！

内容的提问来源于stack exchange，提问作者Ángela Vidal Monge