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

14

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.

Ethan Shepherd
fonte
5
O Weibull é uma distribuição contínua, de modo que a probabilidade de obter exatamente zero tem probabilidade zero. Se você estiver obtendo muitos zeros em seus dados, é uma pista imediata de que o Weibull é inadequado. De qualquer forma, seus dados parecem dados de contagem (ou pelo menos são discretos) e, portanto, um Weibull provavelmente não é a melhor escolha.
cardinal
Adicionar algum contexto sobre a origem dos dados ajudará qualquer um que tente responder tremendamente.
cardinal

Respostas:

8

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


Se os dados acima forem modificados (apenas para ilustração), substituindo os três menores valores por , ou seja,0 0

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)


Edição: Acabei de instalar R para experimentá-lo. Correndo o risco de tornar essa resposta muito longa, para quem estiver interessado, aqui está o meu código R para o método Blischke-Scheuer:

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
res
fonte
4

θfitdistrθθ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θfitdistrθ

jbowman
fonte
2

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.

AdamO
fonte
O OP explicou que um investigador anterior escolheu para caber uma distribuição Weibull, mesmo que os dados são reais contagens - contagens inteiro não negativo do número de ocorrências de algo. Não está claro como seus três casos se relacionam com essa situação.
res
Oh, boa nota! Ajustar-se à distribuição Weibull é flagrantemente errado. Ele tem suporte contínuo e nunca é usado para modelar contagens, mas sim tempos de sobrevivência. Distribuições binomiais negativas seriam uma espécie de distribuição equivalente a dois parâmetros para contagens de modelagem, o que obviamente depende da natureza do processo de geração de dados (do qual temos 0 informações, conforme o problema é declarado). Obrigado por apontar isso para mim.
Adamo
1

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.

John Bauer
fonte
0

[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.

zbicyclist
fonte