Graph of empirical and theoretical distributions for Zero Inflated Poisson Distribution

Question

Following is a kind of data set I am working on it:

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)

I have used the fitdistrplus package to fit this data using Poisson distribution.

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

Once I use the plot() function on x:

plot(x)

Two side-by-side plots are constructed by default. One graph consists of an empirical vs theoretical distribution graph and another has CDF.

The data set has (29/244)=12% zero values thus My concern here is I want to fit the data using Zero Inflated Poisson distribution and want to have the same type of graph that is produced using fitdistrplus package (as I mentioned in earlier part: plot(x))

Answer 1

You could use the estimate of lambda and put it into VGAM::dzipois and 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))

Then make the plot by hand.

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)

Note: There are some other distributions in VGAM that might be worth looking at.