Posbinom: Positive-Binomial Distribution

Description

Density, distribution function, quantile function and random generation for the positive-binomial distribution.

Usage

dposbinom(x, size, prob, log = FALSE)
pposbinom(q, size, prob, lower.tail = TRUE, log.p = FALSE)
qposbinom(p, size, prob, lower.tail = TRUE, log.p = FALSE)
rposbinom(n, size, prob)

Arguments

x, q

vector of quantiles.

vector of probabilities.

number of observations. If length(n) > 1 then the length is taken to be the number required.

size

number of trials. It is the $N$ symbol in the formula given in posbinomial.

prob

probability of success on each trial.

log, log.p, lower.tail

Arguments that are passed on to pbinom etc.

Value

dposbinom gives the density, pposbinom gives the distribution function, qposbinom gives the quantile function, and rposbinom generates random deviates.

Details

The positive-binomial distribution is a binomial distribution but with the probability of a zero being zero. The other probabilities are scaled to add to unity. The mean therefore is $$\mu / (1-(1-\mu)^N)$$ where $\mu$ is the argument prob above. As $\mu$ increases, the positive-binomial and binomial distributions become more similar. Unlike similar functions for the binomial distribution, a zero value of prob is not permitted here.

Examples

Run this code

prob = 0.2; size = 10
y = rposbinom(n = 1000, size, prob)

table(y)
mean(y)  # Sample mean
prob / (1-(1-prob)^size) # Population mean

(ii = dposbinom(0:size, size, prob))
cumsum(ii) - pposbinom(0:size, size, prob)  # Should be 0s
table(rposbinom(100, size, prob))

table(qposbinom(runif(1000), size, prob))
round(dposbinom(1:10, size, prob) * 1000) # Should be similar

barplot(rbind(dposbinom(x = 0:size, size, prob),
                        dbinom(x = 0:size, size, prob)),
        beside = TRUE, col = c("blue","green"),
        main=paste("Positive-binomial(", size, ",", prob, ") (blue) vs",
        " Binomial(", size, ",", prob, ") (green)", sep=""),
        names.arg = as.character(0:size), las=1)

# Simulated data example
nn = 1000; sizeval1 = 10; sizeval2 = 20
pdat <- data.frame(x2 = seq(0, 1, length = nn))
pdat <- transform(pdat, prob1  = logit(-2 + 2 * x2, inverse = TRUE),
                        prob2  = logit(-1 + 1 * x2, inverse = TRUE),
                        sizev1 = rep(sizeval1, len = nn),
                        sizev2 = rep(sizeval2, len = nn))
pdat <- transform(pdat, y1 = rposbinom(nn, size = sizev1, prob = prob1),
                        y2 = rposbinom(nn, size = sizev2, prob = prob2))
with(pdat, table(y1))
with(pdat, table(y2))
# Multivariate response
fit2  = vglm(cbind(y1, y2) ~ x2, posbinomial(mv = TRUE),
             trace  = TRUE, pdat, weight = cbind(sizev1, sizev2))
coef(fit2, matrix = TRUE)