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?