Como ajustar uma distribuição Weibull para inserir dados contendo zeros?

Estou tentando reproduzir um algoritmo de previsão existente, transmitido por um pesquisador aposentado. O primeiro passo é ajustar alguns dados observados a uma distribuição Weibull, para obter uma forma e escala que serão usadas para prever valores futuros. Estou usando R para fazer isso. Aqui está um exemplo do meu 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')

Isso funciona bem, a menos que haja zeros na matriz de entrada, o que faz com que ela falhe completamente. O mesmo acontece no SAS. Pelo que entendi, isso ocorre porque uma das etapas no cálculo da distribuição Weibull é usar o log natural, que é indefinido para 0. Existe uma maneira razoável de contornar isso?

O melhor que encontrei até agora é adicionar 1 a todos os meus valores de entrada, ajustar a curva e subtrair um dos meus valores previstos ("desloque" a curva para cima e depois volte para 1). Isso se encaixa razoavelmente bem nos dados previstos anteriormente, mas parece que deve ser uma maneira errada de fazê-lo.

editar: os valores na matriz de entrada são observados, dados do mundo real (o número de ocorrências de alguma coisa) por um intervalo de anos. Então, em alguns anos, o número de ocorrências foi zero. Seja a melhor maneira ou não (concordo que pode não ser), o autor original do algoritmo afirma ter usado a distribuição Weibull, e eu tenho que tentar replicar o processo deles.

(Como outros salientaram, uma distribuição Weibull provavelmente não é uma aproximação apropriada quando os dados são apenas números inteiros. O objetivo a seguir é apenas ajudá-lo a determinar o que o pesquisador anterior fez, com ou sem razão.)

Existem vários métodos alternativos que não são afetados por zeros nos dados, como o uso de vários estimadores de método de momentos. Isso normalmente requer solução numérica de equações envolvendo a função gama, porque os momentos da distribuição Weibull são dados em termos dessa função. Eu não estou familiarizado com R, mas aqui está um programa Sage que ilustra um dos métodos mais simples - talvez ele possa ser adaptado ao R? (Você pode ler sobre esse e outros métodos em, por exemplo, "A distribuição Weibull: um manual" de Horst Rinne, p. 455ff - no entanto, há um erro de digitação em sua eq.12.4b, como o '-1' é 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)

Isso produziu a saída

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]

então o mesmo procedimento produz a saída

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))
}

Isso reproduz (até cinco dígitos significativos) os dois exemplos do Sábio acima:

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
}

Em seguida, minimize esta função usando a otimização unidimensional:

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

onde acabei de criar o "+10" com base em nada.

Para os dados com os três menores valores substituídos por zeros, obtemos:

> 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$

Deveria falhar, você deveria agradecer que falhou.

Suas observações mostraram que as falhas ocorreram no momento em que você começou a observá-las. Se este é um processo real, proveniente de dados reais (e não de dados simulados), você precisa explicar de alguma forma o motivo pelo qual está obtendo zeros. Eu vi estudos de sobrevivência em que 0 vezes aparecem como consequência de uma das várias coisas:

1. Os dados estão realmente truncados: os objetos estavam em risco e falharam antes do início do estudo e você deseja fingir que os observou o tempo todo.
2. Os instrumentos estão mal calibrados: você não tem precisão de medição suficiente para o estudo e, portanto, as falhas que ocorrem perto da hora de início foram codificadas exatamente como zero.
3. A coisa codificada como zero não é zero. São pessoas ou objetos que foram excluídos da análise de uma maneira ou de outra. O zero apenas aparece nos dados como consequência da mesclagem, classificação ou recodificação de valores ausentes.

Portanto, no caso 1: você precisa usar métodos de censura adequados, mesmo que isso signifique extrair retrospectivamente os registros. O caso 2 significa que você pode usar o algoritmo EM porque possui um problema de precisão. Os métodos bayesianos funcionam de forma semelhante aqui também. O caso 3 significa que você só precisa excluir os valores que deveriam estar ausentes.

Concordo com a resposta do cardeal acima. No entanto, também é bastante comum adicionar uma constante para evitar zeros. Outro valor comumente usado é 0,5, mas qualquer constante positiva pode ter sido usada. Você pode tentar uma variedade de valores para ver se consegue identificar o valor exato usado pelo pesquisador anterior. Então você pode ter certeza de que é capaz de reproduzir os resultados dele antes de procurar uma distribuição melhor.

[Supondo que Weibull seja apropriado] O livro de Johnson Kotz e Balakrishnan tem muitas maneiras de estimar os parâmetros de Weibull. Algumas delas não dependem dos dados que não incluem zeros (por exemplo, usando a média e o desvio padrão, ou certos percentis).

Johnson, NL, Kotz, S. e Balakrishnan, N. (1994). Distribuições univariadas contínuas. Nova York: Wiley, aproximadamente na página 632.