You said you already know $\mu$ and $\sigma$, and you're trying to find $\theta$, so you're solving the equation
$$\sigma^2 = \mu + \mu^2/\theta$$ for $\theta$ with everything else given.
This is elementary algebraic manipulation. I'm not sure I see the difficulty in isolating $\theta$ here.
However, moment-matching is not necessarily a particularly efficient way to estimate parameters in the negative binomial. If you have all the data (not just the mean and standard deviation), why would you not use maximum likelihood to get $\theta$?
This is straightforward in R in several ways. One is to use
library(MASS) #this comes with R but is not loaded by default
summary(glm.nb(y~1))
Here's an example on a set of data:
7 3 7 7 2 6 3 12 12 8 4 6 26 21 8 15 12 13 2 11 0 27
7 8 5 13 2 12 4 1 5 6 17 2 3 6 12 12 2 12
Using the moment-matching approach you get $\hat{\mu}=8.525$ and $\hat{\theta}=2.268387$. The ML estimate of the mean parameter will be the same
but the estimate of $\theta$ is not:
> nbfity = glm.nb(y~1)
> summary(nbfity)
[snip a few lines]
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 2.1430 0.1168 18.35 <2e-16 ***
(Dispersion parameter for Negative Binomial(2.3352) family taken to be 1)
Null deviance: 43.291 on 39 degrees of freedom
Residual deviance: 43.291 on 39 degrees of freedom
AIC: 251.26
Number of Fisher Scoring iterations: 1
Theta: 2.335
Std. Err.: 0.672
2 x log-likelihood: -247.257
So $\hat{\theta}=2.335$ and you get a standard error for free. The estimate of $\mu$ is obtained by exponentiating the coefficient of the intercept:
exp(nbfity$coefficients[1])
[1] 8.525
One can get a standard error on this scale if needed using predict
.
Note that besides MASS::rnegbin
there's stats::rnbinom
which uses another parameterization (which the Wikipedia article on the negative binomial helpfully gives the mean and variance in terms of).