Use a Chi-squared test as explained at https://stats.stackexchange.com/a/17148/919. The R code below implements such a test, with defaults appropriate for calving data.
The Chi-squared test is appropriate for such discrete datasets, as explained at that link. To see that it might work to decide whether a particular dataset is consistent with a Negative Binomial distribution, here are the results of simulating a thousand datasets of 180 independent values.
The first two datasets are shown in the scatterplot at left (the pairings are arbitrary). A histogram of the Chi-squared p-values is shown next. Its small deviations from a uniform distribution (shown by the horizontal gray line) are attributable to chance, strongly indicating this test provides the correct p-values when the null hypothesis (of a Negative Binomial distribution) is true.
The power of this test is its ability to discriminate Negative Binomial from other distributions. For typical calving data, the Negative Binomial is close to Normal (allowing for rounding to the nearest day). So are other distributions, such as Poisson distributions with appropriate parameters. Thus, we shouldn't expect much of this test (or any such test). The distributions of p-values from simulated data with Poisson and Normal distributions appear in the right two histograms. Because there is a tendency for p-values to be smaller with these alternatives, the test has some power to detect the difference. But because the p-values aren't very small very often, the power is low: with a dataset of 180, it will be difficult to distinguish Negative Binomial from Poisson from Normal data. This suggests that the question whether the data are consistent with a Negative Binomial distribution might have little inherent meaning or usefulness.
The parameters for this example come from Werth, Azzam and Kinder, Calving intervals in beef cows at 2, 3, and 4 years of age when breeding is not restricted after calving. J. Animal Sci. 1996, 74:593-596. Because this paper does not provide adequate descriptive statistics, I estimated the mean and variance (and set the breaks for the chi-squared test) from this figure:
This is the R code to implement the calculations and plots shown here. It's not bulletproof: before applying any of these functions to other datasets, it would be prudent to test them and perhaps include code to verify the maximum likelihood estimates are correct.
library(MASS) # rnegbin
#
# Specify parameters to generate data.
#
mu <- 360 # Mean days in interval
v <- 30^2 # Variance of days: must exceed mu^2
n <- 18000 # Sample size
n.sim <- 3e2 # Simulation size
#
# Functions to fit a negative binomial to data.
#
pnegbin <- function(k, mu, theta) {
v <- mu + mu^2/theta # Variance
p <- 1 - mu / v # "Success probability"
r <- mu * (1-p) / p # "Number of failures until the experiment is stopped"
pbeta(p, k+1, r, lower.tail=FALSE)
}
# #
# # Test `pnegbin` by comparing it to randomly generated data.
# #
# z <- rnegbin(1e3, mu, theta)
# plot(ecdf(z))
# curve(pnegbin(x, mu, theta), add=TRUE, col="Red", lwd=2)
#
# Maximum likelihood fitting of data based on counts in predefined bins.
# Returns the fit and chi-squared statistics.
#
negbin.fit <- function(x, breaks) {
if (missing(breaks))
breaks <- c(-1, seq(-40, 30, by=10) + 365, Inf)
observed <- table(cut(x, breaks))
n <- length(x)
counts.expected <- function(n, mu, theta)
n * diff(pnegbin(breaks, mu, theta))
log.lik.m <- function(parameters) {
mu <- parameters[1]
theta <- parameters[2]
-sum(observed * log(diff(pnegbin(breaks, mu, theta))))
}
v <- var(x)
m <- mean(x)
if (v > m) theta <- m^2 / (v - m) else theta <- 1e6 * m^2
parameters <- c(m, theta)
fit <- optim(parameters, log.lik.m)
expected <- counts.expected(n, fit$par[1], fit$par[2])
chi.square <- sum(res <- (observed - expected)^2 / expected)
df <- length(observed) - length(parameters) - 1
p.value <- pchisq(chi.square, df, lower.tail=FALSE)
return(list(fit=fit, chi.square=chi.square, df=df, p.value=p.value,
data=x, breaks=breaks, observed=observed, expected=expected,
residuals=res))
}
#
# Test on randomly generated data.
#
# set.seed(17)
sim <- replicate(n.sim, negbin.fit(rnegbin(n, mu, theta))$p.value)
#
# Generate data for illustration.
#
theta <- mu^2 / (v - mu)
x <- rnegbin(n, mu, theta)
y <- rnegbin(n, mu-4.3, theta)
#
# Display data and simulation.
#
par(mfrow=c(1,4))
plot(x-365, y-365, pch=15, col="#00000040",
xlab="First calving interval", ylab="Second calving interval",
main="Simulated Data")
abline(h=0)
abline(v=0)
hist(sim, freq=FALSE, xlab="p-values", ylab="Frequency",
main="Histogram of Simulated P-values",
sub="Negative Binomial Data")
abline(h=1, col="Gray", lty=3)
#
# Simulate non-Negative Binomial data for comparison.
#
sim.2 <- replicate(n.sim, negbin.fit(rpois(n, mu))$p.value)
hist(sim.2, freq=FALSE, xlab="p-values", ylab="Frequency",
main="Histogram of Simulated P-values",
sub="Poisson Data")
abline(h=1, col="Gray", lty=3)
sim.3 <- replicate(n.sim, negbin.fit(floor(rnorm(n, mu, sqrt(mu))))$p.value)
hist(sim.3, freq=FALSE, xlab="p-values", ylab="Frequency",
main="Histogram of Simulated P-values",
sub="Normal Data")
abline(h=1, col="Gray", lty=3)
par(mfrow=c(1,1))
Best Answer
sum(dnbinom(x=0:10,size=10,prob=0.6))is unlikely to be correct as it does not use the $5$ in the question.The help page says "This represents the number of failures which occur in a sequence of Bernoulli trials before a target number of successes is reached" and "
sizetarget for number of successful trials", so to get $5$ successes in up to to $10$ trials, you can afford $0$ through to $5$ failures. So either of these will work (the second uses the cumulative distribution function):You could reformulate the question using a binomial distribution: you want at least $5$ successes in $10$ attempts. So any of these work too: