HelpersMG (version 5.1)

dSnbinom: Density for the sum of random variable with negative binomial distributions.

Description

Density for the sum of random variable with negative binomial distributions. If all prob values are the same, infinite is automatically set to 0.

Usage

dSnbinom(
  x = stop("You must provide a x value"),
  size = NULL,
  prob = NULL,
  mu = NULL,
  log = FALSE,
  tol = 1e-06,
  infinite = 1000
)

Arguments

x

vector of (non-negative integer) quantiles.

size

target for number of successful trials, or dispersion parameter (the shape parameter of the gamma mixing distribution). Must be strictly positive, need not be integer.

prob

probability of success in each trial. 0 < prob <= 1.

mu

alternative parametrization via mean.

log

logical; if TRUE, probabilities p are given as log(p).

tol

Tolerance for recurrence

infinite

Maximum level of recursion

Value

dSnbinom gives the density

Details

dSnbinom returns the density for the sum of random variable with negative binomial distributions

References

Furman, E., 2007. On the convolution of the negative binomial random variables. Statistics & Probability Letters 77, 169-172.

See Also

Other Distribution of sum of random variable with negative binomial distributions: pSnbinom(), qSnbinom(), rSnbinom()

Examples

# NOT RUN {
library(HelpersMG)
alpha <- c(1, 2, 5, 1, 2)
p <- c(0.1, 0.12, 0.13, 0.14, 0.14)
dSnbinom(20, size=alpha, prob=p)
dSnbinom(20, size=alpha, prob=p, log=TRUE)
dSnbinom(20, size=2, mu=c(0.01, 0.02, 0.03))
dSnbinom(20, size=2, mu=c(0.01, 0.02, 0.03), log=TRUE)
# Test with a single distribution
dSnbinom(20, size=1, mu=20)
# when only one distribution is available, it is the same as dnbinom()
dnbinom(20, size=1, mu=20)
# If a parameter is supplied as only one value, it is supposed to be constant
dSnbinom(20, size=1, mu=c(14, 15, 10))
# The function is vectorized:
plot(0:200, dSnbinom(0:200, size=alpha, prob=p), bty="n", type="h", xlab="x", ylab="Density")
# Comparison with simulated distribution using rep replicates
alpha <- c(2.1, 2.05, 2)
mu <- c(10, 30, 20)
rep <- 100000
distEmpirique <- rSnbinom(rep, size=alpha, mu=mu)
tabledistEmpirique <- rep(0, 301)
names(tabledistEmpirique) <- as.character(0:300)
tabledistEmpirique[names(table(distEmpirique))] <- table(distEmpirique)/rep

plot(0:300, dSnbinom(0:300, size=alpha, mu=mu), type="h", bty="n", 
   xlab="x", ylab="Density", ylim=c(0,0.02))
plot_add(0:300, tabledistEmpirique, type="l", col="red")
legend(x=200, y=0.02, legend=c("Empirical", "Theoretical"), 
   text.col=c("red", "black"), bty="n")

# Example with the approximation mu=mean(mu)
plot(0:300, dSnbinom(0:300, size=alpha, mu=mu), type="h", bty="n", 
   xlab="x", ylab="Density", ylim=c(0,0.02))
plot_add(0:300, tabledistEmpirique, type="l", col="red")
legend(x=200, y=0.02, legend=c("Empirical", "Theoretical"), 
   text.col=c("red", "black"), bty="n")
   
# example to fit the distribution
data <- rnbinom(1000, size=1, mu=10)
hist(data)
ag <- rep(1:100, 10)
r <- aggregate(data, by=list(ag), FUN=sum)
hist(r[,2])

parx <- c(size=1, mu=10)

dSnbinomx <- function(x, par) {
  -sum(dSnbinom(x=x[,2], mu=rep(par["mu"], 10), size=par["size"], log=TRUE))
}

fit_mu_size <- optim(par = parx, fn=dSnbinomx, x=r, method="BFGS", control=c(trace=TRUE))
fit_mu_size$par
# }