Ошибка в FitDist с гамма-распределением

Question

Ниже приведены мои коды:

library(fitdistrplus)

s <- c(11, 4, 2, 9, 3, 1, 2, 2, 3, 2, 2, 5, 8,3, 15, 3, 9, 22, 0, 4, 10, 1, 9, 10, 11,
 2, 8, 2, 6, 0, 15, 0 , 2, 11, 0, 6, 3, 5, 0, 7, 6, 0, 7, 1, 0, 6, 4, 1, 3, 5,
 2, 6, 0, 10, 6, 4, 1, 17, 0, 1, 0, 6, 6, 1, 5, 4, 8, 0, 1, 1, 5, 15, 14, 8, 1,
 3, 2, 9, 4, 4, 1, 2, 18, 0, 0, 10, 5, 0, 5, 0, 1, 2, 0, 5, 1, 1, 2, 3, 7)

o <- fitdist(s, "gamma", method = "mle")
summary(o)
plot(o)

---

и ошибка говорит:

Ошибка в fitdist (s, "gamma", method = "mle"): функция mle не смогла оценить параметры, с кодом ошибки 100

Answer 1

Гамма-распределение не допускает нулевых значений (вероятность будет равна нулю, а логарифмическая вероятность будет бесконечной, для ответа 0), если параметр формы не равен точно 1,0 (то есть экспоненциальное распределение - см. ниже) ... это статистическая / математическая проблема, а не проблема программирования. Вам нужно будет найти что-нибудь разумное, чтобы сделать с нулевыми значениями. В зависимости от того, что имеет смысл для вашего приложения, вы можете (например)

• выберите другой дистрибутив для тестирования
• исключить нулевые значения (fitdist(s[s>0], ...))
• установить нулевые значения на некоторое разумное ненулевое значение (fitdist(replace(s,which(s==0),0.1),...)

, который (если таковой имеется) лучше всего зависит от вашей заявки .

---

@ Первый ответ Сандипана Дея (оставляя нули в наборе данных) выглядит логично , но на самом деле он застревает при параметре формы, равном 1.

o <- fitdist(s, "exp", method = "mle")

дает тот же ответ, что и код @ Сандипана (за исключением того, что он оценивает скорость = 0,2161572, обратный к параметру масштаба = 4,626262, который оценивается для гамма-распределения - это просто изменение параметризации). Если вы решили использовать экспоненту вместо гаммы, это нормально, но вы должны делать это нарочно, а не случайно ...

Чтобы проиллюстрировать, что подгонка с нулями может работать не так, как ожидалось, я создам свою собственную функцию отрицательного логарифмического правдоподобия и выведу поверхность вероятности для каждого случая.

mfun <- function(sh,sc,dd=s) {
  -sum(dgamma(dd,shape=sh,scale=sc,log=TRUE))
}

library(emdbook)  ## for curve3d() helper function

Поверхность с нулями:

cc1 <- curve3d(mfun(x,y),
      ## set up "shape" limits" so we evaluate
      ## exactly shape=1.000 ...
        xlim=c(0.55,3.55),
        n=c(41,41),
        ylim=c(2,5),
        sys3d="none")


png("gammazero1.png")
with(cc1,image(x,y,z))
dev.off()

В этом случае поверхность определяется только при форме = 1 (то есть экспоненциальное распределение); белые области представляют бесконечные логарифмические вероятности. Дело не в том, что форма = 1 является лучшей посадкой, а в том, что это только посадка ...

Поверхность, исключенная из нулей:

cc2 <- curve3d(mfun(x,y,dd=s[s>0]),
               ## set up "shape" limits" so we evaluate
               ## exactly shape=1.000 ...
               xlim=c(0.55,3.55),
               n=c(41,41),
               ylim=c(2,5),
               sys3d="none")


png("gammazero2.png")
with(cc2,image(x,y,z))
with(cc2,contour(x,y,z,add=TRUE))
abline(v=1.0,lwd=2,lty=2)
dev.off()

Answer 2

Просто укажите начальные значения параметров гамма-распределения (масштаб, форма), которые нужно вычислить с помощью mle с использованием optim, а также нижние границы для параметров, это должно работать.

o <- fitdist(s, "gamma", lower=c(0,0), start=list(scale=1,shape=1))
summary(o)

#Fitting of the distribution ' gamma ' by maximum likelihood 
#Parameters : 
#      estimate Std. Error
#scale 4.626262         NA
#shape 1.000000         NA
#Loglikelihood:  -250.6432   AIC:  505.2864   BIC:  510.4766 

Согласно комментариям @Ben Bolker, мы можем сначала исключить нулевые точки:

o <- fitdist(s[s!=0], "gamma", method = "mle", lower=c(0,0), start=list(scale=1,shape=1))
summary(o)
#Fitting of the distribution ' gamma ' by maximum likelihood 
#Parameters : 
#      estimate Std. Error
#scale 3.401208         NA
#shape 1.622378         NA
#Loglikelihood:  -219.6761   AIC:  443.3523   BIC:  448.19