График эмпирических и теоретических распределений для нулевого завышенного распределения Пуассона

Question

Ниже приводится набор данных, над которым я работаю:

      data <- c(0, 1, 0, 11, 2, 0, 3, 0, 0, 2, 1, 3, 1, 0, 1, 0, 0, 0, 2, 3, 
0, 0, 0, 8, 1, 1, 1, 0, 1, 1, 2, 7, 0, 0, 0, 5, 2, 3, 6, 1, 1, 
5, 2, 9, 0, 0, 1, 21, 16, 2, 9, 6, 25, 2, 1, 12, 16, 14, 15, 
15, 6, 1, 12, 12, 13, 5, 5, 6, 4, 7, 11, 8, 4, 5, 8, 3, 8, 4, 
7, 4, 7, 2, 5, 6, 4, 5, 1, 0, 8, 5, 6, 8, 9, 8, 9, 7, 7, 9, 8, 
9, 4, 4, 7, 13, 9, 13, 12, 10, 9, 8, 7, 11, 5, 5, 0, 1, 33, 4, 
22, 19, 22, 9, 5, 4, 17, 7, 7, 4, 5, 3, 0, 0, 9, 3, 0, 0, 36, 
40, 5, 4, 0, 11, 0, 7, 5, 25, 39, 26, 4, 20, 12, 4, 17, 3, 22, 
12, 14, 8, 9, 11, 7, 11, 10, 9, 16, 6, 24, 8, 5, 6, 14, 3, 9, 
4, 1, 20, 0, 1, 7, 9, 0, 12, 2, 29, 56, 16, 8, 28, 0, 19, 25, 
35, 87, 56, 66, 60, 58, 14, 10, 12, 13, 13, 34, 26, 18, 13, 22, 
13, 12, 15, 41, 11, 11, 11, 5, 6, 7, 8, 8, 17, 16, 12, 21, 38, 
34, 10, 77, 41, 7, 12, 1, 16, 20, 8, 5, 2, 20, 7, 16, 12, 6, 
10, 31, 12)

Я использовал fitdistrplus пакет для соответствия этим данным с использованием распределения Пуассона.

      library(fitdistrplus)
x <- fitdist(data, "pois")

Как только я использую plot() функционировать на x:

      plot(x)

По умолчанию строятся два расположенных рядом графика. Один график состоит из графика эмпирического и теоретического распределения, а другой - из CDF.

Набор данных имеет (29/244)=12% нулевых значений, поэтому меня беспокоит то, что я хочу подогнать данные с использованием распределения Пуассона с нулевым надуванием и хочу иметь тот же тип графика, который создается с использованием пакета fitdistrplus (как я упоминал в ранняя часть: plot(x))

Answer 1

Вы можете использовать оценку лямбды и поместить ее в VGAM::dzipois а также VGAM::pzipois.

      library(VGAM)
distr1 <- dzipois(0:max(data), lambda=x$estimate, pstr0=mean(data == 0))
distr2 <- pzipois(0:max(data), lambda=x$estimate, pstr0=mean(data == 0))

Затем сделайте сюжет своими руками.

      op <- par(mfrow=c(1, 2))
plot(proportions(table(data)), ylim=c(0, max(proportions(table(data)))), 
     main="Empirical and \ntheoretical distribution", ylab='Density')
lines((0:max(data)) + 0.2, distr1, type="h", lwd=2, col=2)
legend("topright", c("empirical", "theoretical"), lwd=2, col=1:2, cex=0.8)
plot(ecdf(data), pch=1, col.hor=0, main="Empirical and \ntheoretial CDFs")
lines(0:max(data), distr2, type="s", lwd=2, col=2)
legend("bottomright", c("empirical", "theoretical"), lwd=2, col=1:2, cex=0.8)
par(op)

Примечание. В VGAM есть и другие дистрибутивы, на которые стоит обратить внимание.