¿Cómo ajustar una distribución de Weibull para ingresar datos que contienen ceros?

Estoy tratando de reproducir un algoritmo de predicción existente, transmitido por un investigador retirado. El primer paso es ajustar algunos datos observados a una distribución de Weibull, para obtener una forma y una escala que se utilizarán para predecir valores futuros. Estoy usando R para hacer esto. Aquí hay un ejemplo de mi código:

x<-c(23,19,37,38,40,36,172,48,113,90,54,104,90,54,157,51,77,78,144,34,29,45,16,15,37,218,170,44,121)
f<-fitdistr(x, 'weibull')

Esto funciona bien a menos que haya ceros en la matriz de entrada, lo que hace que falle por completo. Lo mismo sucede en SAS. Según tengo entendido, esto se debe a que uno de los pasos para calcular la distribución de Weibull es tomar el registro natural, que no está definido para 0. ¿Hay alguna forma razonable de evitar esto?

Lo mejor que he encontrado hasta ahora es sumar 1 a todos mis valores de entrada, ajustar la curva y luego restar uno de mis valores predichos ("desplazar" la curva hacia arriba y luego hacia abajo por 1). Esto se ajusta bastante bien a los datos previstos anteriormente, pero parece que debe ser una forma incorrecta de hacerlo.

editar: los valores en la matriz de entrada se observan, datos del mundo real (el número de ocurrencias de algo) durante un rango de años. Entonces, en algunos años, el número de ocurrencias fue cero. Ya sea que sea la mejor manera o no (estoy de acuerdo en que puede no serlo), el autor del algoritmo original afirma haber usado la distribución Weibull, y tengo que intentar replicar su proceso.

(Como otros han señalado, es probable que una distribución de Weibull no sea una aproximación apropiada cuando los datos son solo enteros. Lo siguiente tiene la intención de ayudarlo a determinar lo que hizo el investigador anterior, correcta o incorrectamente).

Existen varios métodos alternativos que no se ven afectados por ceros en los datos, como el uso de varios estimadores de métodos de momentos. Por lo general, requieren una solución numérica de ecuaciones que involucran la función gamma, porque los momentos de la distribución de Weibull se dan en términos de esta función. No estoy familiarizado con R, pero aquí hay un programa Sage que ilustra uno de los métodos más simples: ¿tal vez se pueda adaptar a R? (Puede leer sobre este y otros métodos similares en, por ejemplo, "La distribución de Weibull: un manual" de Horst Rinne, p. 455ff; sin embargo, hay un error tipográfico en su ecuación 12.4b, como el '-1' es redundante)

"""
Blischke-Scheuer method-of-moments estimation of (a,b)
for the Weibull distribution F(t) = 1 - exp(-(t/a)^b)
""" 

x = [23,19,37,38,40,36,172,48,113,90,54,104,90,54,157,
      51,77,78,144,34,29,45,16,15,37,218,170,44,121]
xbar = mean(x)
varx = variance(x)
var("b"); f(b) = gamma(1+2/b)/gamma(1+1/b)^2 - 1 - varx/xbar^2
bhat = find_root(f, 0.01, 100)
ahat = xbar/gamma(1+1/bhat)
print "Estimates: (ahat, bhat) = ", (ahat, bhat)

Esto produjo la salida

Estimates: (ahat, bhat) =  (81.316784310814455, 1.3811394719075942)

x = [23,0,37,38,40,36,172,48,113,90,54,104,90,54,157,
      51,77,78,144,34,29,45,0,0,37,218,170,44,121]

entonces el mismo procedimiento produce la salida

Estimates: (ahat, bhat) =  (78.479354097488923, 1.2938352346035282)

fit_weibull <- function(x)
{
    xbar <- mean(x)
    varx <- var(x)
    f <- function(b){return(gamma(1+2/b)/gamma(1+1/b)^2 - 1 - varx/xbar^2)}
    bhat <- uniroot(f,c(0.02,50))$root
    ahat <- xbar/gamma(1+1/bhat)
    return(c(ahat,bhat))
}

Esto reproduce (a cinco dígitos significativos) los dos ejemplos de Sage anteriores:

x <- c(23,19,37,38,40,36,172,48,113,90,54,104,90,54,157,
     51,77,78,144,34,29,45,16,15,37,218,170,44,121)
fit_weibull(x)
[1] 81.316840  1.381145

x <- c(23,0,37,38,40,36,172,48,113,90,54,104,90,54,157,
      51,77,78,144,34,29,45,0,0,37,218,170,44,121)
fit_weibull(x)
[1] 78.479180  1.293821

$\theta$fitdistr$\theta$$\theta$fitdistr

foo <- function(theta, x)
{
  if (theta <= -min(x)) return(Inf);
  f <- fitdistr(x+theta, 'weibull')
  -2*f$loglik
}

Luego minimice esta función utilizando la optimización unidimensional:

bar <- optimize(foo, lower=-min(x)+0.001, upper=-min(x)+10, x=x)

donde acabo de inventar el "+10" basado en nada en absoluto.

Para los datos con los tres valores más pequeños reemplazados por ceros, obtenemos:

> bar
$minimum
[1] 2.878442

$objective
[1] 306.2792

> fitdistr(x+bar$minimum, 'weibull')
     shape        scale   
   1.2836432   81.1678283 
 ( 0.1918654) (12.3101211)
> 

bar$minimum$\theta$fitdistr$\theta$

Debería fallar, deberías estar agradecido de que haya fallado.

Sus observaciones mostraron que las fallas ocurrieron en el mismo momento en que comenzó a observarlas. Si este es un proceso real, proveniente de datos reales (y no de datos simulados), debe de alguna manera explicar la razón por la que obtiene ceros. He visto estudios de supervivencia donde aparecen 0 veces como consecuencia de una de varias cosas:

1. Los datos están realmente truncados: los objetos estaban en riesgo y fallaron antes de que comenzara el estudio y usted quiere fingir que los ha observado todo el tiempo.
2. Los instrumentos están mal calibrados: no tiene suficiente precisión de medición para el estudio y, por lo tanto, las fallas que ocurren cerca de la hora de inicio se codificaron exactamente como cero.
3. La cosa codificada como cero no es un cero. Son personas u objetos que fueron excluidos del análisis de una forma u otra. El cero solo aparece en los datos como consecuencia de la fusión, clasificación o recodificación de los valores faltantes.

Entonces, para el caso 1: debe usar métodos de censura adecuados, incluso si eso significa extraer registros retrospectivamente. El caso 2 significa que puede usar el algoritmo EM porque tiene un problema de precisión. Los métodos bayesianos también funcionan de manera similar aquí. El caso 3 significa que solo necesita excluir los valores que se suponía que faltaban.

Estoy de acuerdo con la respuesta del cardenal anterior. Sin embargo, también es bastante común agregar una constante para evitar ceros. Otro valor comúnmente usado es 0.5, pero podría haberse usado cualquier constante positiva. Puede probar un rango de valores para ver si puede identificar el valor exacto utilizado por el investigador anterior. Entonces podría estar seguro de que puede reproducir sus resultados, antes de buscar una mejor distribución.

[Asumiendo que Weibull es apropiado] El libro de Johnson Kotz y Balakrishnan tiene muchas maneras de estimar los parámetros de Weibull. Algunos de estos no dependen de los datos que no incluyen ceros (por ejemplo, usando la media y la desviación estándar, o usando ciertos percentiles).

Johnson, NL, Kotz, S. y Balakrishnan, N. (1994). Distribuciones Univariadas Continuas. Nueva York: Wiley, aproximadamente en la página 632.