Parameterization of negative binomial densities (dnbinom) in R

Question

I hope it will be here, not a math forum.

I am trying to calculate the negative binomial density in R. My parameters are currently in the form mu(average number of expected failures) and the parameter is overdispersion k. Since I call dnbinomin the compiled C code, I need to convert these parameters to the default parameterization dnbinom, which expects the probability of success in each test pand the number of successes nuntil the count stops, I used the following equations to solve pboth nin terms of muand k:

n = mu*p/(1-p)
mu + mu^2/k = n(1-p)/p^2

After a bit of algebra I get

p = k/mu + 1
n = -(mu^2+k*mu)/k

Checking these equations with various parameterizations dnbinomin R reveals the problem:

> k = 1.2
> mu = 15
> p = k/mu+1
> n = -(mu*k+mu^2)/k
> dnbinom(10,size=n,prob=p)
[1] NaN
Warning message:
In dnbinom(x, size, prob, log) : NaNs produced
> dnbinom(10,mu=mu,size=k)
[1] 0.03560668

What about the parameterization of R that I will skip? I am pretty sure my algebra is correct. Alternatively, is there a way to stick to the original parameterization (in terms of muand k) when called dnbinomfrom C?

Answer 1

As they ?dnbinomsay: the size " must be strictly positive , it is not necessary to be whole"

In your example

 n = -(mu*k+mu^2)/k
 n
 [1] -202.5

So that dnbinomdoes not work.

I think that n should be mu*p/(1-p)after the calculation p(as you wrote above), therefore as follows:

k = 1.2
mu = 15
p = k/(k+mu)
n = mu*p/(1-p)
dnbinom(10,size=n,prob=p)
    [1] 0.03560668
dnbinom(10,mu=mu,size=k)
    [1] 0.03560668

Which looks good.