Jak dopasować rozkład Weibulla do danych wejściowych zawierających zera?

Próbuję odtworzyć istniejący algorytm prognozowania, przekazany przez emerytowanego badacza. Pierwszym krokiem jest dopasowanie niektórych obserwowanych danych do rozkładu Weibulla, aby uzyskać kształt i skalę, które zostaną wykorzystane do przewidywania przyszłych wartości. Używam do tego R. Oto przykład mojego kodu:

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

Działa to dobrze, chyba że w tablicy wejściowej znajdują się zera, co powoduje całkowite niepowodzenie. To samo dzieje się w SAS. Jak rozumiem, dzieje się tak, ponieważ jednym z kroków w obliczaniu rozkładu Weibulla jest pobranie logarytmu naturalnego, który jest niezdefiniowany dla 0. Czy istnieje rozsądny sposób na obejście tego?

Najlepsze, co do tej pory znalazłem, to dodać 1 do wszystkich moich wartości wejściowych, dopasować krzywą, a następnie odjąć jedną z moich przewidywanych wartości („przesunąć” krzywą w górę, a następnie w dół o 1). To dość dobrze pasuje do wcześniej przewidywanych danych, ale wydaje się, że musi to być zły sposób.

edycja: Wartości w tablicy wejściowej są obserwowane w rzeczywistych danych (liczba wystąpień czegoś) przez szereg lat. Tak więc w niektórych latach liczba wystąpień wynosiła zero. Niezależnie od tego, czy jest to najlepszy sposób, czy nie (zgadzam się, że tak nie jest), autor oryginalnego algorytmu twierdzi, że użył dystrybucji Weibulla i muszę spróbować powtórzyć ich proces.

(Jak zauważyli inni, rozkład Weibulla prawdopodobnie nie będzie odpowiednim przybliżeniem, gdy dane są tylko liczbami całkowitymi. Poniższe informacje mają na celu pomóc ci ustalić, co zrobił poprzedni badacz, słusznie lub niesłusznie).

Istnieje kilka alternatywnych metod, na które zer danych nie ma wpływu, takich jak stosowanie różnych metod estymatorów momentów. Zazwyczaj wymagają one numerycznego rozwiązania równań obejmujących funkcję gamma, ponieważ momenty rozkładu Weibulla podane są w kategoriach tej funkcji. Nie znam R, ale oto program Sage ilustrujący jedną z prostszych metod - może można go dostosować do R? (O tym i innych takich metodach można przeczytać w, np. „Dystrybucja Weibulla: podręcznik” Horsta Rinne'a, s. 455ff - jednak w jego równaniu 12.4b jest literówka, jako „-1” jest zbędny).

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

To dało wynik

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]

wtedy ta sama procedura daje wynik

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

To odtwarza (do pięciu cyfr znaczących) dwa powyższe przykłady Sage:

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
}

Następnie zminimalizuj tę funkcję, stosując optymalizację jednowymiarową:

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

gdzie właśnie stworzyłem „+10” w oparciu o nic.

W przypadku danych z trzema najmniejszymi wartościami zastąpionymi zerami otrzymujemy:

> 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 jest MLE z $\theta$, a dane fitdistrwyjściowe to MLE parametrów Weibulla, łącznie z$\theta$to jest. Jak widać, są one dość zbliżone do estymatorów metody momentów @res pokazanych powyżej.

Powinien zawieść, powinieneś być wdzięczny, że zawiódł.

Twoje obserwacje wykazały, że awarie miały miejsce w momencie, gdy zacząłeś je obserwować. Jeśli jest to prawdziwy proces pochodzący z rzeczywistych (a nie symulowanych danych), musisz jakoś wyjaśnić, dlaczego otrzymujesz zera. Widziałem badania dotyczące przeżycia, w których 0 razy pojawia się jako konsekwencja jednej z kilku rzeczy:

1. Dane są w rzeczywistości obcinane: obiekty były zagrożone i zawiodły przed rozpoczęciem badania i chcesz udawać, że je obserwowałeś przez cały czas.
2. Instrumenty są źle skalibrowane: nie masz wystarczającej dokładności pomiaru do badania, więc awarie występujące w pobliżu czasu rozpoczęcia zostały zakodowane jako dokładnie zero.
3. Rzecz zakodowana jako zero nie jest zerem. Są to osoby lub przedmioty, które zostały wykluczone z analizy w ten czy inny sposób. Zero pojawia się tylko w danych w wyniku scalenia, sortowania lub w inny sposób przekodowania brakujących wartości.

Tak więc w przypadku 1: musisz zastosować odpowiednie metody cenzury, nawet jeśli oznacza to retrospektywne pobieranie zapisów. Przypadek 2 oznacza, że ​​możesz użyć algorytmu EM, ponieważ masz problem z precyzją. Metody bayesowskie również tutaj działają podobnie. Przypadek 3 oznacza, że ​​musisz po prostu wykluczyć wartości, których brakowało.

Zgadzam się z powyższą odpowiedzią kardynała. Jednak dość często dodaje się stałą, aby uniknąć zer. Inną powszechnie stosowaną wartością jest 0,5, ale można było zastosować dowolną stałą dodatnią. Możesz wypróbować szereg wartości, aby sprawdzić, czy możesz dokładnie określić wartość użytą przez poprzedniego badacza. Wtedy możesz mieć pewność, że jesteś w stanie odtworzyć jego wyniki, zanim zaczniesz szukać lepszej dystrybucji.

[Zakładając, że Weibull jest odpowiedni] Książka Johnsona Kotza i Balakrishnana ma wiele sposobów na oszacowanie parametrów Weibulla. Niektóre z nich nie zależą od danych nie zawierających zer (np. Przy użyciu średniej i odchylenia standardowego lub przy użyciu określonych percentyli).

Johnson, NL, Kotz, S. i Balakrishnan, N. (1994). Ciągłe rozkłady jednowymiarowe. Nowy Jork: Wiley, mniej więcej na stronie 632.