Estimating shared parameter from negative binomial sample 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?
 A: 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.
A: 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.
