parameterization of the negative binomial in scipy via mean and std

Question

I am trying to fit my data to a Negative Binomial Distribution with the package scipy in Python. However, my validation seems to fail.

These are my steps:

1. I have some demand data which is described by the statistics:

mu = 1.4
std = 1.59
print(mu, std)

2. I use the parameterization function below, taken from this post to compute the two NB parameters.

def convert_params(mu, theta):
    """
    Convert mean/dispersion parameterization of a negative binomial to the ones scipy supports

    See https://en.wikipedia.org/wiki/Negative_binomial_distribution#Alternative_formulations
    """
    r = theta
    var = mu + 1 / r * mu ** 2
    p = (var - mu) / var
    return r, 1 - p

I pass (hopefully correctly...) my two statistics - the naming convention between different sources is rather confusing at this point p, r, k

firstParam, secondParam = convert_params(mu, std)

3. I would then use these two parameters to fit the distribution:

from scipy.stats import nbinom

rv = nbinom(firstParam, secondParam)

Then I calculate a value R with the Percent Point Function .ppf(0.95). The value R in the context of my problem is a Reorder Point.

R = rv.ppf(0.95)

4. Now is when I expect to validate the previous steps, but I do not manage to retrieve my original statistics mu and std with mean and math.sqrt(var) respectively.

import math

mean, var = nbinom.stats(firstParam, secondParam, moments='mv')
print(mean, math.sqrt(var))

What am I missing? Any feedback about the parameterization implemented in Scipy?

Immortelle

Answer 1

Conversion code is wrong, I believe, SciPy is NOT using Wiki convention, but Mathematica convention

#%%
import numpy as np
from scipy.stats import nbinom

def convert_params(mean, std):
    """
    Convert mean/dispersion parameterization of a negative binomial to the ones scipy supports

    See https://mathworld.wolfram.com/NegativeBinomialDistribution.html
    """
    p = mean/std**2
    n = mean*p/(1.0 - p)
    return n, p

mean = 1.4
std  = 1.59

n, p = convert_params(mean, std)

print((n, p))

#%%

m, v = nbinom.stats(n, p, moments='mv')
print(m, np.sqrt(v))

Code prints back 1.4, 1.59 pair

And reorder point computed as

rv = nbinom(n, p)
print("reorder point:", rv.ppf(0.95))

outputs 5

Greatuncle

Answer 2

It looks like you are using a different conversion. The last bullet at the cited wikipedia section gives the formulas shown below. With these formulas you get back the exact same mu and std:

import numpy as np
from scipy.stats import nbinom

def convert_mu_std_to_r_p(mu, std):
    r = mu ** 2 / (std ** 2 - mu)
    p = 1 - mu / std ** 2
    return r, 1 - p

mu = 1.4
std = 1.59
print("mu, std:", mu, std)
firstParam, secondParam = convert_mu_std_to_r_p(mu, std)
mean, var = nbinom.stats(firstParam, secondParam, moments='mv')
print("mean, sqrt(var):", mean, np.sqrt(var))

rv = nbinom(firstParam, secondParam)
print("reorder point:", rv.ppf(0.95))

Output:

mu, std: 1.4 1.59
mean, sqrt(var): 1.4 1.59
reorder point: 5.0

Acetum