Trudności z dopasowaniem rozkładu gamma za pomocą R

Question

Próbuję oszacować parametry dopasowania rozkładu gamma do gęstości ekologicznej (tj. Biomasy na obszar). Używam polecenia fitdistr() z pakietu MASS w R (wersja 3.0.0: x86_64-w64-mingw32/x64 (64-bit)). Jest to polecenie szacowania maksymalnego prawdopodobieństwa dla parametrów dystrybucji.

Wektory danych są dość duże, ale statystyki podsumowujące są następujące:

Min. = 0; 1st Qu. = 87.67; Median = 199.5; Mean = 1255; Variance = 2.79E+07; 3rd Qu. = 385.6; Max. = 33880

Kod używam, aby uruchomić procedurę MLE jest:

gdist <- fitdistr(data, dgamma, 
        start=list(shape=1, scale=1/(mean(data))),lower=c(1,0.1)) 

R daje mi się następujący błąd:

Error in optim(x = c(6.46791148085828, 4060.54750836902, 99.6201565968665, : non-finite finite-difference value [1]

Inni, którzy doświadczyli tego typu problemów i zwrócił się do Stackov erflow na pomoc najwyraźniej znalazł rozwiązanie dodając argument "lower =" do ich kodu i/lub usuwając zera. Uważam, że R dostarczy parametry dla dopasowania, jeśli usunę obserwacje zerowe, ale miałem wrażenie, że rozkłady gamma obejmowały zakres 0 < = x> inf (Forbes i wsp. 2011. Rozkład statystyczny)?

Czy mam złe zdanie na temat zasięgu dystrybucji gamma? Czy jest jakiś inny problem, którego mi brakuje odnośnie MLE (w którym jestem nowicjuszem).

Answer 1

Pierwsze szacunkowe metodą momentów (pasująca się średnią = kształt * skalę i wariancji = kształt * skalę^2) mamy

mean <- 1255 
var <- 2.79e7 
shape = mean^2/var ## 0.056 
scale = var/mean  ## 22231 

Teraz wygenerować pewne dane z tej dystrybucji:

set.seed(101) 
x = rgamma(1e4,shape,scale=scale) 
summary(x) 
##  Min. 1st Qu. Median  Mean 3rd Qu.  Max. 
##  0.00  0.00  0.06 1258.00  98.26 110600.00 

MASS::fitdistr(x,"gamma") ## error 

Podejrzewam, że problem polega na tym, że wywołanie optim zakłada, że ​​parametry (kształt i skala lub kształt i szybkość) są w przybliżeniu tej samej wielkości, które nie są. Można obejść ten problem przez skalowanie danych:

(m <- MASS::fitdistr(x/2e4,"gamma")) ## works fine 
##  shape   rate  
## 0.0570282411 0.9067274280 
## (0.0005855527) (0.0390939393) 

fitdistr daje parametr szybkości zamiast parametru skali: aby wrócić do parametru kształtu chcesz odwrócić i ponownie skalę ...

1/coef(m)["rate"]*2e4 ## 22057 

Nawiasem mówiąc, fakt, że kwantyle danych symulowanych nie pasują do twoich danych bardzo dobrze (np. Mediana z x = 0,06 wobec mediany 199 w danych) sugerują, że twoje dane mogą nie pasować do Gamma wszystkie to dobrze - np może być kilka wartości odstających wpływających na średnią i wariancję bardziej niż kwantyle?

PS powyżej użyłem wbudowany „gamma” prognozy w fitdistr zamiast używać dgamma: z wartości począwszy oparty na metodzie momentów i skalowanie danych przez 2E4, to działa (choć daje ostrzeżenie o NaNs produced chyba że określić lower)

m2 <- MASS::fitdistr(x/2e4,dgamma, 
     start=list(shape=0.05,scale=1), lower=c(0.01,0.01))