Repository: scipy/scipy
Pull request: DOC: stats: note another common parameterization of nbinom
#### Reference issue
gh-11465 (does not close)
Closes gh-13287
Checks a box in gh-4538

#### What does this implement/fix?
SciPy parameterized the negative binomial distribution in terms of the  probability of success `p`. Another common parameterization of the negative binomial distribution is in terms of the mean number of failures before achieving `n` successes. This has come up in the referenced issues, so I have added the conversion information to the documentation.

Also, I noticed that the docs were building with `n=0.4` as the number of successes. Apparently, `nbinom` accepts non-integer values of `n`, and `n=0.4` was intentionally added as a test case in `_distn_params.py`. I'm going to leave that alone, but the documentation should show an integer value of `n`. I've rearranged the two test cases to hopefully build the docs with `n=5`.

Problem statement:
Generate binomial CDF with mu instead of prob
**Is your feature request related to a problem? Please describe.**

The feature is related to a problem where I am unclear as to how to translate my data into parameters: `n` and `p` for `nbinom.cdf`. I am sort of lost translating R to Python here given different params:

```
test_arr = np.array([(1, 1, 3, 1, 1, 0, 0, 2, 2, 0, 1, 3, 0, 0, 0, 0, 3, 2, 1, 5, 2, 2,
           4, 0, 7, 3, 0, 0, 2, 1, 2, 0, 1, 4, 1, 0, 0, 0, 0, 0, 0, 0, 1, 1,
           0, 0, 0, 0, 0, 2, 0, 3, 3, 1, 3, 0, 0, 2, 0, 3, 1, 0, 0, 0, 0, 1,
           0, 4, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 2, 1, 2, 0, 0, 3, 2, 2, 0, 3,
           0, 1, 0, 1, 0, 0, 0, 2, 1, 0, 1, 0)])
nbinom.cdf(test_arr, n=?, p=?)  # how determine n and p from the data?
```




**Describe the solution you'd like**

When using R, I can obtain `nbinom.cdf` seemingly from the data itself
```
library(pscl)
arr1 <- c(1, 1, 3, 1, 1, 0, 0, 2, 2, 0, 1, 3, 0, 0, 0, 0, 3, 2, 1, 5, 2, 2,
           4, 0, 7, 3, 0, 0, 2, 1, 2, 0, 1, 4, 1, 0, 0, 0, 0, 0, 0, 0, 1, 1,
           0, 0, 0, 0, 0, 2, 0, 3, 3, 1, 3, 0, 0, 2, 0, 3, 1, 0, 0, 0, 0, 1,
           0, 4, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 2, 1, 2, 0, 0, 3, 2, 2, 0, 3,
           0, 1, 0, 1, 0, 0, 0, 2, 1, 0, 1, 0)
m1 <- zeroinfl(arr1 ~ 1, dist='negbin')
m1$theta  # dispersion param

library(VGAM)
# size is dispersion parameter
pnbinom(arr1, size=m1$theta, mu=mean(arr1))
```

It would be nice to derive the CDF of a value/vector by specifying the dispersion parameter and the mean.


**Describe alternatives you've considered**
 N/A

**Additional context (e.g. screenshots)**
 N/A

Hints from issue discussion:
statsmodels is supposed to have some functions to convert from GLM/Regression parameterization of negbin to the standard negbin.

scipy avoids duplicating distributions with just changes in parameterization
Wow thanks again for all your time and help. Are you able to help me better understand how I can determine `n` and `p`, if say, my data was the numpy array: `test_arr `?
you can look at the code of NegativeBinomialP at predict and convert_params methods, for appropriate version of "P"

I don't find convert_params for NegativeBinomial right now, but old issue and recipe is in
https://github.com/statsmodels/statsmodels/issues/106#issuecomment-43961704
@afogarty85 When you have figured this out, would you let us know what we should add to the documentation to help others? Just recently I noticed a similar question in gh-11465. 
Thanks all, this has been super helpful, particularly this: https://github.com/statsmodels/statsmodels/issues/106#issuecomment-43961704.

In terms of documentation, perhaps it might be helpful to add @josef-pkt's work showing how we can estimate `n` and `p` if we already have collected some data and then we want to proceed to other methods like the CDF.

```
from scipy.stats import nbinom
import statsmodels.api as sm
import numpy as np

nbinom_rv = nbinom.rvs(n=1.089, p=0.514, size=10000)
intercept = np.ones(len(nbinom_rv))
res = sm.NegativeBinomial(nbinom_rv, intercept, loglike_method='nb2').fit();
mu = np.exp(res.params[0])
alpha = res.params[1]
res.params

# find size and prob
Q = 0  # nb2, cited most common
size = 1. / alpha * mu ** Q
prob = size / (size + mu)
print('estimated params:', 'size:', round(size, 4), 'prob:', round(prob, 4))
cdf_out = nbinom.cdf(nbinom_rv, n=size, p=prob)
np.allclose(nbinom_rv, nbinom.ppf(cdf_out, size, prob))
```







@afogarty Looking more closely, part of your question is about fitting the distribution to data. That is a feature SciPy doesn't have for discrete distributions, but it has already been requested in gh-11948. In the meantime, you have statsmodels and, if you're interested, a [recipe](https://github.com/scipy/scipy/issues/11948#issuecomment-742964862) in gh-11948. Since fitting a  discrete distribution to data is not part of SciPy right now, though, we don't mention it in the documentation for any discrete distribution, and so I don't think we should add information about that just to `nbinom`.

The other part is about parameterization in terms of mu, and I think it is worth including that information in the documentation. Would gh-13292 address that and close this PR? Please post your thoughts there. Thanks!

Commit message: DOC: stats: fixup in nbinom

Patch to review:
diff --git a/scipy/stats/_discrete_distns.py b/scipy/stats/_discrete_distns.py
index 1fe26de158e8..b76e9d57a5f0 100644
--- a/scipy/stats/_discrete_distns.py
+++ b/scipy/stats/_discrete_distns.py
@@ -31,8 +31,8 @@ class binom_gen(rv_discrete):
     for ``k`` in ``{0, 1,..., n}``, :math:`0 \leq p \leq 1`
 
     `binom` takes ``n`` and ``p`` as shape parameters,
-    where :math:`p` is the probability of a single success 
-    and :math:`1-p` is the probability of a single failure. 
+    where :math:`p` is the probability of a single success
+    and :math:`1-p` is the probability of a single failure.
 
     %(after_notes)s
 
@@ -113,8 +113,8 @@ class bernoulli_gen(binom_gen):
     for :math:`k` in :math:`\{0, 1\}`, :math:`0 \leq p \leq 1`
 
     `bernoulli` takes :math:`p` as shape parameter,
-    where :math:`p` is the probability of a single success 
-    and :math:`1-p` is the probability of a single failure. 
+    where :math:`p` is the probability of a single success
+    and :math:`1-p` is the probability of a single failure.
 
     %(after_notes)s
 
@@ -258,8 +258,29 @@ class nbinom_gen(rv_discrete):
     for :math:`k \ge 0`, :math:`0 < p \leq 1`
 
     `nbinom` takes :math:`n` and :math:`p` as shape parameters where n is the
-    number of successes, :math:`p` is the probability of a single success, 
-    and :math:`1-p` is the probability of a single failure. 
+    number of successes, :math:`p` is the probability of a single success,
+    and :math:`1-p` is the probability of a single failure.
+
+    Another common parameterization of the negative binomial distribution is
+    in terms of the mean number of failures :math:`\mu` to achieve :math:`n`
+    successes. The mean :math:`\mu` is related to the probability of success
+    as
+
+    .. math::
+
+       p = \frac{n}{n + \mu}
+
+    The number of successes :math:`n` may also be specified in terms of a
+    "dispersion", "heterogeneity", or "aggregation" parameter :math:`\alpha`,
+    which relates the mean :math:`\mu` to the variance :math:`\sigma^2`,
+    e.g. :math:`\sigma^2 = \mu + \alpha \mu^2`. Regardless of the convention
+    used for :math:`\alpha`,
+
+    .. math::
+       \begin{eqnarray*}
+       p & = & \frac{\mu}{\sigma^2} \\
+       n & = & \frac{\mu^2}{\sigma^2 - \mu}
+       \end{eqnarray*}
 
     %(after_notes)s
 
@@ -327,8 +348,8 @@ class geom_gen(rv_discrete):
 
     for :math:`k \ge 1`, :math:`0 < p \leq 1`
 
-    `geom` takes :math:`p` as shape parameter, 
-    where :math:`p` is the probability of a single success 
+    `geom` takes :math:`p` as shape parameter,
+    where :math:`p` is the probability of a single success
     and :math:`1-p` is the probability of a single failure.
 
     %(after_notes)s
@@ -1106,7 +1127,7 @@ def _rvs(self, a, size=None, random_state=None):
         #     https://www.sciencedirect.com/science/
         #     article/abs/pii/S0378375804003519
         # Furthermore, the two-sided geometric distribution is
-        # equivalent to the difference between two iid geometric 
+        # equivalent to the difference between two iid geometric
         # distributions.
         #   Reference (page 179):
         #     https://pdfs.semanticscholar.org/61b3/
diff --git a/scipy/stats/_distr_params.py b/scipy/stats/_distr_params.py
index c04c25560d20..78b69edf3448 100644
--- a/scipy/stats/_distr_params.py
+++ b/scipy/stats/_distr_params.py
@@ -122,8 +122,8 @@
     ['hypergeom',(21,3,12)],  # numpy.random (3,18,12) numpy ticket:921
     ['hypergeom',(21,18,11)],  # numpy.random (18,3,11) numpy ticket:921
     ['logser', (0.6,)],  # re-enabled, numpy ticket:921
-    ['nbinom', (5, 0.5)],
     ['nbinom', (0.4, 0.4)],  # from tickets: 583
+    ['nbinom', (5, 0.5)],
     ['planck', (0.51,)],   # 4.1
     ['poisson', (0.6,)],
     ['randint', (7, 31)],


Code review comments:
diff_hunk: @@ -258,8 +258,29 @@ class nbinom_gen(rv_discrete):
     for :math:`k \ge 0`, :math:`0 < p \leq 1`
 
     `nbinom` takes :math:`n` and :math:`p` as shape parameters where n is the
-    number of successes, :math:`p` is the probability of a single success, 
-    and :math:`1-p` is the probability of a single failure. 
+    number of successes, :math:`p` is the probability of a single success,
+    and :math:`1-p` is the probability of a single failure.
+
+    Another common parameterization of the negative binomial distribution is
+    in terms of the mean number of failures :math:`\mu` to achieve :math:`n`
+    successes. The mean :math:`\mu` is related to the probability of success
+    as
+
+    .. math::
+
+       p = \frac{n}{n + \mu}
+
+    The number of successes :math:`n` may also be specified in terms of a
+    "dispersion", "heterogeneity", or "aggregation" parameter :math:`\alpha`,
+    which relates the mean :math:`\mu` to the variance :math:`\sigma^2`,
+    e.g. :math:`\sigma^2 = \mu + \alpha \mu^2`. Regardless of the convention
+    used for :math:`\alpha`,
+
+    .. math::
+       \begin{eqnarray*}
+       p & = & \frac{\mu}{\sigma^2} \\
+       n & = & \frac{\mu^2}{\sigma^2 - \mu}
+       \end{eqnarray*}
path: scipy/stats/_discrete_distns.py
text: @author:
```suggestion

       p &= \frac{\mu}{\sigma^2} \\
       n &= \frac{\mu^2}{\sigma^2 - \mu}

```