Estimating shared parameter from negative binomial sample

Question

I came across a problem and I think it is unsolvable, but I would like to make sure that this is the case.

Let consider a sample: $$ X_1, X_2, ... X_n \sim NB(r_i, p). $$ Therefore I know that observations share $p$ parameter, but each observation $X_i$ can have a different $r_i$. I do not know anything about $r_i$. Is it possible to estimate the common parameter $p$? I do not need to estimate $r_i$.

I tried to find maximum likelihood estimator (forgetting that $r_i$ differs through the sample) and retrieve $p$, but it gives wrong result.

Here I assume that $p=0.2$.

r <- rgamma(1000, 3,0.1)
x <- sapply(r, function(x) rnbinom(1, size=x, p=0.2))

library(MASS)
par.nb<- fitdistr(x, "negative binomial")
par.nb$estimate

size <- par.nb$estimate[1]
mu <- par.nb$estimate[2]

p_est <-size/(mu+size)
p_est

But $p\_est = 0.02085524$

Any suggestions?

Answer 1

fitdistr finds you an estimate for $r$ and $p$ as if all the data is sampled from the same distribution. I you try to find an MLE for $p$ assuming all $r_i$ are different, you'll have $$\hat p = \frac{\sum x_i}{\sum (x_i+r_i)}$$

However, you also don't know $r_i$. If you try to take derivative wrt $r_i$, you'll encounter something like $$\psi(x_i+r_i)-\psi(r_i)+\log\left({r_i \over {r_i+x_i}}\right)=0$$ where $\psi(x)$ is digamma function. As wikipedia page says, there is no closed form solution exist. Moreover, for estimating $r_i$, you're using only one sample. But, as cited in the wiki page:

The maximum likelihood estimator only exists for samples for which the sample variance is larger than the sample mean.

When $n=1$, the sample variance is $0$, so we can't actually estimate $r_i$. So, a MLE is not possible.

Answer 2

It is obviously impossible to estimate $r_1,\ldots,r_n$ and $p$ because that is $n+1$ parameters and you have only $n$ data points from which to estimate them. It is impossible to estimate $n+1$ independent parameters from $n$ data values.

You can however estimate $p$ if the $r_i$ are known or if you at least know $\bar r=\frac1n\sum_{i=1}^n r_i$. In the latter case, the sample mean $\bar x$ would be an unbiased estimator of $$\frac{\bar r p}{1-p}$$ so $$\hat p=\frac{\exp(\bar x/r)}{1+\exp(\bar x/r)}$$ is a consistent estimator of $p$.

This estimator $\hat p$ is somewhat inefficient however if the $r_i$ are very different. If you did know the values of $r_i$ then a more efficient estimator of $p$ would be possible.