munb <- 5; size <- 4; n <- 1000
table(y <- rposnegbin(n, munb=munb, size=size))
mean(y) # sample mean
munb / (1 - (size/(size+munb))^size) # population mean
munb / pnbinom(0, mu=munb, size=size, lower.tail = FALSE) # same as before
x <- (-1):17
(ii <- dposnegbin(x, munb=munb, size=size))
max(abs(cumsum(ii) - pposnegbin(x, munb=munb, size=size))) # Should be 0
x <- 0:10
barplot(rbind(dposnegbin(x, munb=munb, size=size), dnbinom(x, mu=munb, size=size)),
beside = TRUE, col = c("blue","green"),
main=paste("dposnegbin(munb=", munb, ", size=", size, ") (blue) vs",
" dnbinom(mu=", munb, ", size=", size, ") (green)", sep=""),
names.arg = as.character(x))
# Another test for pposnegbin()
nn <- 5000
mytab <- cumsum(table(rposnegbin(nn, munb=munb, size=size))) / nn
myans <- pposnegbin(sort(as.numeric(names(mytab))), munb=munb, size=size)
max(abs(mytab - myans)) # Should be 0
Run the code above in your browser using DataLab