# Find confidence interval for negative binomial using Mathematica

I would like to find the exact 95% confidence intervals on the mean of a negative binomial random variable using Mathematica. It is possible to approximate the distribution using a gamma/normal, but not for parameter ranges in which I am interested. As such, I would like a way of numerically deriving the 95% confidence intervals

The way in which Mathematica defines the negative binomial distribution, is in terms of $n$ and $p$, and the mean can be calculated as:

Mean[NegativeBinomialDistribution[n, p]]


which returns $\frac{n(1-p)}{p}$. It is this estimated quantity which I would like to construct confidence intervals around.

Does anyone have an idea as to how to start here? I must admit I am a bit stuck as to where to start.

Best,

Ben

• You need to be a bit clearer about your purpose. Do you need the interval within which 95% of the m variates drawn from NegativeBinomialDistribution[n, p] fall? Or do you want the CI around the measured mean describing the 95% probability of the position of their real mean? – Sjoerd C. de Vries Mar 25 '15 at 13:57
• I want a 95% confidence interval around the estimated mean. Does that make sense? I have edited the question now to reflect this. Best, Ben – ben18785 Mar 25 '15 at 13:59

The distribution of the means of sufficiently large sets of random drawings from a given distribution (well, many of them) approaches a normal distribution with a mean equal to the mean of the given distribution and a standard deviation equal to the standard deviation of the given distribution divided by the square root of the sample size -1.

Using this information, you obtain the theoretical bounds using Quantile:

Block[{sampleSize = 100, p = 0.1, lb = 0.025, ub = 0.975, mean},
mean = Mean[NegativeBinomialDistribution[sampleSize, p]];
Quantile[
NormalDistribution[
mean,
StandardDeviation[NegativeBinomialDistribution[sampleSize, p]]/
Sqrt[sampleSize - 1]
],
{lb, ub}
]
]

(* {881.312, 918.688} *)


A simulation:

Block[{sampleSize = 100, p = 0.1, simSampleSize = 1000000, lb = 0.025, ub = 0.975},
Quantile[
Mean /@
RandomVariate[
NegativeBinomialDistribution[sampleSize,p],
{simSampleSize, sampleSize}
]//N,
{lb, ub}
]
]
(* {881.47, 918.7} *)


Update: Exact approach

For low values of the sample size the normal distribution approximation doesn't hold. One can obtain the numerical exact (but unfortunately not parametrizable) result as follows:

Block[{sampleSize = 5, n = 9, p = 0.1, lb = 0.025, ub = 0.975},
Quantile[
TransformedDistribution[
Total[Table[x[i], {i, sampleSize}]],
Table[x[i], {i, sampleSize}] \[Distributed]
NegativeBinomialDistribution[n, p]
]
], {lb, ub}
]/sampleSize
] // N

(* {57.8, 107.6} *)


Note the division by sampleSize after the Quantile function. This is done because Quantile can handle the Total of the random variables, but not the Mean or Total[...]/sampleSize. The result of ranking the means should be the same as the ranking of totals followed by a division by a constant.

Compare with simulation:

Block[{sampleSize = 5, n = 9, p = 0.1, simSampleSize = 100000, lb = 0.025, ub = 0.975},
Quantile[
Mean /@
RandomVariate[
NegativeBinomialDistribution[n, p],
{simSampleSize, sampleSize}
] // N,
{lb, ub}]
]

(* {57.8, 107.4} *)


How it works:

What I'm doing is to define a distribution that is the sum of sampleSize stochastic variables each distributed as a NegativeBinomialDistribution. This mimics the sample process that you perform if you draw from this distribution. This transformed distribution too has a PDF, CDF, InverseCDF and because of the latter also a Quantile function. Mathematica is able to calculate this function and so I can find the 2.5 and 97.5 percentile levels which define the 95% confidence interval. As I explained above, I used the Total with a division later on instead of the Mean directly as Mathematica doesn't seem to handle that in this case.

Necessity of exact approach

To show under which conditions the normal distribution approach would be sufficient we can generate a few plots using various levels of n, p and sample size. Note the erratic behavior of the exact predictions that match the simulations very well. Not so for the normal distribution approach for low sample sizes. The rule of thumb that one can use the approximation for sample sizes of 30 and up is confirmed.

ParallelTable[
With[{
l1 = If[n == 64,
PlotLegends -> PointLegend[{"Normal approach"}], {}],
l2 = If[n == 64, PlotLegends -> PointLegend[{"Exact"}], {}],
l3 = If[n == 64, PlotLegends -> PointLegend[{"Simulation"}], {}]},
sims =
Table[Prepend[simulate[sampleSize, n, p, 100000, 0.025, 0.975],
sampleSize], {sampleSize, {3, 5, 8, 13, 21, 34}}];
Show[
DiscretePlot[
normalApp[sampleSize, n, p, 100000, 0.025, 0.975], {sampleSize,
3, 35}, PlotStyle -> {Dashed, Blue}, Filling -> None,
Joined -> True, l1],
DiscretePlot[
exactApp[sampleSize, n, p, 100000, 0.025, 0.975], {sampleSize, 3,
35}, PlotStyle -> {Dashed, Darker@Green}, Filling -> None,
Joined -> True, l2],
ListPlot[{sims[[All, {1, 2}]], sims[[All, {1, 3}]]},
PlotStyle -> Red, l3]
, Frame -> True, Axes -> None,
FrameLabel -> {If[p == 0.95, "Sample size", None],
If[n == 2, "2.5 / 97.5 percentile", None],
StringTemplate["n = . p ="][n, p], ""}, ImageSize -> 250
]
], {p, {0.1, 0.5, 0.8, 0.95}}, {n, {2, 8, 64}}
] // Grid


• Just take out the Quantile block from the first block of code and evaluate that. It will give you the analytic result. – Sjoerd C. de Vries Mar 25 '15 at 14:35
• @ben see update – Sjoerd C. de Vries Mar 25 '15 at 15:00
• @ben18785 tiny update: separated sample size and the parameter n. – Sjoerd C. de Vries Mar 25 '15 at 17:25
• @ben Sorry, no paper refs, but it's actually not that complex. I added a bit of explanation above. – Sjoerd C. de Vries Mar 25 '15 at 17:49
• @ben18785 And one final update... – Sjoerd C. de Vries Mar 25 '15 at 22:58