Трудность подбора гамма-распределения с R

Я пытаюсь оценить параметры для гамма-распределения, соответствующие данным экологической плотности (т.е. по биомассе на область). Я использовал команду fitdistr() из пакета MASS в R (версия 3.0.0: x86_64-w64-mingw32/x64 (64-разрядная версия)). Это команда оценки максимального правдоподобия для параметров распределения.

Векторы данных довольно велики, но сводная статистика выглядит следующим образом:

Мин. = 0; 1-й кв. = 87,67; Медиана = 199,5; Среднее value = 1255; Разница = 2.79E + 07; 3 кв. = 385,6; Максимум. = 33880

Код, который я использую для запуска процедуры MLE:

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

R дает мне следующую ошибку:

Ошибка в оптимизации (x = c (6.46791148085828, 4060.54750836902, 99.6201565968665,: конечное конечное разностное значение [1]

Другие, которые испытали этот тип проблемы и обратились к stackoverflow для справки, похоже, нашли решение в добавлении аргумента "lower =" к своему коду и/или удалению нулей. Я обнаружил, что R будет предоставлять параметры для подгонки, если я удалю нулевые наблюдения, но у меня создалось впечатление, что гамма-распределения охватывают диапазон 0 <= x > inf (Forbes et al., 2011. Статистические распределения)?

Я получил неправильное впечатление относительно диапазона гамма-распределения? Или есть еще одна проблема, которую я пропускаю в отношении MLE (в которой я новичок).

Ответы

Ответ 1

Получение приблизительной оценки методом моментов (совпадение среднего масштаба = форма * шкала и дисперсия = форма * шкала ^ 2) мы имеем

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

Теперь создайте некоторые данные из этого дистрибутива:

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

Я сильно подозреваю, что проблема заключается в том, что базовый вызов optim предполагает, что параметры (форма и масштаб, форма и скорость) имеют примерно такую ​​же величину, какими они не являются. Вы можете обойти это, масштабируя свои данные:

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

fitdistr дает параметр скорости, а не параметр масштаба: чтобы вернуться к требуемому параметру формы, инвертировать и повторно масштабировать...

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

Кстати, тот факт, что квантили симулированных данных не очень хорошо соответствуют вашим данным (например, медиана x= 0,06 против медианы 199 в ваших данных), предположили, что ваши данные могут не соответствовать Гамма все это хорошо - например могут быть некоторые выбросы, влияющие на среднее и дисперсию больше, чем квантили?

PS выше Я использовал встроенную оценку "гамма" в fitdistr вместо использования dgamma: со стартовыми значениями, основанными на методе моментов, и масштабированием данных на 2e4, это работает (хотя это дает предупреждение о NaNs produced, если мы не укажем lower)

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